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ABSTRACT 

A threshold detection system based on the Neyman-Pearson 
criterion is derived for an infrared nutating optical system. 
Detection is optimum for small signal- to-noise ratios and a 
Gaussian uncertainty in the pointing error of the optical 
system. The signal spectrum and background power spectral 
density for a rectangular space filter are computed numeri- 
cally and used to specify the matched filter for the thresh- 
old detection system and evaluate its performance. 
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I. INTRODUCTION 



An infrared (IR) detection system determines the 
presence or absence of a target in a particular area of 
space by processing the received IR energy from that area. 

The target must be found amid ever-present background 

/ 

radiance and any internal noise generated by the detection 
system. In this thesis, one type of system, called a 
nutating optical system, is considered. 

Infrared detection is a combination of spatial frequency 
processing and temporal processing. Spatial frequency process- 
ing is a means of achieving a modulated output from the IR 
radiance passing into the optical system so that information 
may be easily extracted [Ref. 14]. To extract this informa- 
tion, some form of temporal processing is used. In a 
nutating optical system, the optical axis rotates about an 
axis perpendicular to the image plane in such a way that a 
point image traces out a circle in the image plane. The 
spatial frequency filter is a piece of IR sensitive material 
in the image plane which produces a voltage output that is 
amplified and fed into a temporal filter for processing. 

A nutating optical system is frequently used in a missile 
seeker because it produces a signal when there is no 
tracking error and the rotating optics can be made into a 
gyro that provides the inertial reference. 

This work is a synthesis and extension of two recently 
published articles in the literature. The first, by 
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Harger [Ref. 3] , formulates the detection theory of a 
known signal in background and white noise. Since in IR 
target detection the signal is never known completely, 
Harger' s work is extended to the case where the target 
amplitude and position are unknown. However, a known 
target shape is assumed as well as a priori statistics of 
amplitude and position. Since detection is most difficult 
when the amplitude is small, a test (called the threshold 
detector) is derived that is optimum for weak signals. 

To specify the form of the threshold detector and to 
calculate its performance, an eigenvalue integral equation 
must be solved. Ordinarily, this is a difficult undertaking 
but the solution becomes ' trivial if the covariance function 
of the background process is periodic. This is exactly 
the condition realized for the background process out of 
a nutating detector provided that the phase, which is 
physically unimportant, is averaged out. 

Samuelsson's results [Ref. 11] specify the form of the 
signal and background coefficients in terms of the system 
parameters and signal and background radiance. These 
derivations were first published in a Swedish internal 
report [Ref. 10], then re-derived in English by Yoshitani 
[Ref. 6]. The latter's derivations are included in the 
appendices . 

Samuelsson's equations have been applied to a rectangu- 
lar IR detector assuming a Wiener spectrum for the random 
background and the coefficients have been calculated. 
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Using these coefficients, the optimum filter has been 
determined and the probability of detection computed as 
a function of signal-to-noise ratios and background- to- 
noise ratios. 
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II. SPATIAL PROCESSING 



The spatial frequency filtering portion of the detection 
system is comprised of a lens system and an IR sensitive 
detector located in the image plane, the focal plane 
of the optics. To adequately derive information 

from the IR radiation, the lens system and detector must 
be of a nature that causes the voltage output to be modulated 
in a way that information about a target within the field 
of view can be processed. In a practical nutating optical 
system, the lens system moves and the detector is stationary. 
However, for conceptual purposes, this is equivalent to a 
stationary lens systen and nutating detector. 



A. TARGET SPECTRUM 

The object to be detected, the target, is located 
in the object plane a distance R from the lens and parallel 
to the image plane as shown in Figure 1.1. The angle from 
the perpendicular axis connecting the planes to the target 
coordinates measured in the x-direction is called x and is 
measured in milliradians ; the angle in the y-direction is 
y. The target is considered to be close to the optical 
axis and at a sufficient distance, R, such that small angle 
approximations may be made, 



x = tan x 
y = tan y 



( 2 . 1 ) 
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Figure 1.1. Optical Imagin 



Consider a target with a radiance distribution 
2 

N(r) (W/cm -sr) which is the radiant power into a unit solid 
angle per unit projected area of the sources. r =r(x,y) is 
a vector whose origin is the perpendicular axis connecting 
the coordinates . 

Assuming no atmospheric attenuation (this can be 
included in N(r) if it exists) and perfect optics, the power 
distribution in the image plane is 

Pi(r) = N(r) • A 0 (w/sr) (2.2) 

where Aq is the effective entrance area of the optics. 

Dividing by Aq , one can see that the radiance 
distribution in the image, plane is numerically equal to the 
radiance distribution in the object plane. 

However, because the optics is not perfect, a 
point object given by 

Hp 5 (r) = Hp 5 (x) 5 (y) (w) , (2. ) 

2 

(where Hp (W/cm ) is the irradiance at the optics received 
from the point source) is imaged as a power distribution 

Hp • -f (r) (W/cm^ Sr) (2.4) 

where fp(r) (sr ^) is the point spread function of the 
optics. Irradiance is defined as the radiant flux indicent 
on a surface of unit area [Ref. 5]. 

Therefore, the radiance distribution, N(r), in the 
object plane will be imaged as 

N'(r) = J N(r-s) £q(s) d^s (W/cm^Sr) (2.5) 
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2 

where d s = dxdy and the integration is over the entire 
image plane. 

The detector is a space filter in the image plane 
which may be a moving retical or simply an aperture across 
which the power is integrated. 

The power in the image plane incident on the detec- 
tor is 

H = Jn ' ( r) x (r) d 2 r (W/cm 2 ) (2.6) 

where x(r) = x(x,y) is the transmittance of the detector. 

When the detector moves as a function of time, the 
power incident on the detector becomes a function of time, 

H(t) = Jn'(t) x(r-p(t)) d 2 r (W/cm 2 ) (2.7) 

where p(t) describes the movement of the detector in the 
image plane. 

A nutating detector moves circularly in the image 
plane but its orientation remains fixed, i.e., the coordinate 
system (x',y') shown in Fig. 1.2 does not rotate. The nuta- 
tion radius measured with respect to the optical axis is 
given by p. 

Assuming a nutating optics and no relative motion 
between the target and the optics, the radiance on the 
detector will be periodic and thus can be expanded in a 
Fourier series 

00 j nw t 

H ( t ) = 2 H e 0 0 < t < T (2.8) 

n 



13 



where 



"n ‘ f J H(t) 



-J n 0) o t 

e dt . 



Replacing H(t) by equation (2.7), 



, rT -jnw t 

H =^j dte ° d z r N' (r) (r-p(t)) 



(2.9) 



o 



The motion describing function, p(t), for circular 
nutation is 



p(t) = p ( c o s o) Q t , sin w o t) , 



( 2 . 10 ) 



where w Q = 2 tt/T (rad/sec) is the radian frequency of nutation 
and T is the period for one nutation. 

Substituting (2.10) in (2.9), can be shown (see 
Appendix A) to be 



H n = j n /n ’ (k) T *00 e' jn<f) J n (27rp Vk x 2 + k y 2 ) d 2 k 



( 2 . 11 ) 
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Figure 1.2. Circular Nutation in the Image Plane. 



where k = (k , k ) , two-dimensional spacial Fourier 
~ x y 

transform coordinates [Ref. 8], 

N ' (k) = transform of radiance distribution in the 

image plane, 

x* (k) = complex conjugate of the transform of the 

optical transmittance function in the detector 
coordinate system, 

<{> = tan 1 

p = nutation radius, 

and J (z) is the nth order Bessel function of the first 
n 

kind. It has been assumed that the image plane is very 
large compared to the radiance distribution from a target 
so that the image plane may be considered infinite for 
mathematical purposes. 

B. CORRELATION OF BACKGROUND NOISE 

The background power incident on the detector is from 
a sample background scene where a scene is a two-dimensional 
random process characterized by the radius vector r = ( x >y) 

from the optical axis. 

Let N(r) be a sample scene radiance distribution on the 

object plane from an ensemble of scenes which has been given a 

suitable probability structure. The background power incident 
on the detector is 

B(t) = / N' (r) x(r-p(t)) d r . (2.12) 

The covariance function is defined by taking the expected 
value 



IS 



(2.13) 



E{ [ B (t) } - E{ B (t) } ] [B(u) - E{ B (u) } ) = R fi (t,u) 

where E denotes the expectation over the ensemble of scenes. 
Because the detector is nutating, the random process, B(t), 
will not be stationary. However, if the detector and 
background do not move with respect to each other, then 
R B (t^,t 2 ) is doubly periodic: 

R B (t 1 + MT, t 2 + NT) = R fi (t 1 ,t 2 ) (2.14) 



where M and N are integers. 

It can be shown [Ref. 7] that the form of the covari- 
ance of a doubly period process is 



R B ( t ; 






oo 

z 



lmw t, 
J o 1 



mn 



jnoj 0 t 2 



(2.15) 



Because the nutation is periodic, over the ensemble of 
scenes, there is little significance in the starting time t. 
Therefore, considering t as a random variable that is 
uniform over the nutation period, one can "average it out" 
of the covariance function as 

1 T 

R b (t) = i / Rg (t+x , t) dt . (2.16) 

Substituting (2.16) in (2.15), 

00 jmw . T j (m-n) to t 

■ 2 * m e ° T { e ° dt - ' 2 - 17 ^ 

m ,n 0 
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Using the relation 



i T j (u-v)co n t 

6uv — 77 ? f Q dt , (2.18) 

1 0 

the correlation function becomes 

°° jnco n T 

R b (t) = Z 3 n e . (2.19) 

n 

Equation (2,19) shows that B(t) can be considered a 
wide sense stationary process and that it is periodic in 
the mean square sense; i.e., 

R^(x) = Rg(T+x) for any t (2.20) 

The power spectral density is the temporal Fourier 
transform of R^(t). Taking the transform one has 

oo 

Sr>( t0 ) = 2 tt E 8 <S(co-no 3 n ). ( 2 . 21 ) 

B n n 0 

The coefficients, 3 , can be related to the optical 

* n r 

system parameters by 

B n = /^(k)| 2 |F 0 (k)| 2 W B (k) J n 2 (2up / k x 2 +k y 2 ) d 2 k (2.22) 

where (k) is the transform of the point spread function, 
and Wg(k) is the Wiener spectrum, the transform of the back- 
ground correlation function. This derivation was first made 
by Samuelsson [Ref. 10] in an internal Swedish report then 
again by Yoshitani [Ref. 6]. Yoshanti's derivation has 
been included as Appendix B. 

] 7 



Equation (2.19) implies that the background, B(t), can 
be represented by a Fourier series with uncorrelated 
coefficients [Ref. 7]. Thus, 



B (t) = E b e 
n n 



j n u)q t 



(2.23) 



where the convergence is in the mean square sense. 
The coefficients, b , satisfy 



E(b n ) = 



E(B(t)) 

0 



n = 0 
n f 0 



(2.24) 



Moreover, they are uncorrelated 



E (b b ) = 
v in n 



3 



n 



m = n 



0 m f n 



(2.25) 



C. WHITE NOISE COMPONENT 

In an actual system there will be a certain amount of 
internal noise generated in the detector and preamplifier 
which is additive to any background noise or signal. In 
the frequencies of interest the primary source of this 
noise is thermal agitation. This type of noise, called 
Johnson noise, is assumed to be white and Gaussian. 



Let W(t) be Gaussian white noise with power spectral 
N 2 

density y (V /H ). The power spectrum at the output of the 
^ z 

preamp with bandwidth B is shown in Figure 1.3. 
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Figure 1.3. Power Spectrum of White Noise. 

The RMS voltage measured at the output of the preamp is 



V 2 rms = | • 2B = NB (V 2 ) . 



(2.26) 
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III. TEMPORAL PROCESSING 



The output voltage from the detector is temporally pro- 
cessed to determine the presence of a target. The type of 
processor is determined by applying statistical detection 
theory to the known behavior of the nutating system and var- 
ious realistic assumptions about the system and the target. 

An optimum processor in the sense of small signal amplitudes, 
called a threshold detector, is developed in this section. 

Harger [Ref. 3] derived the detection theory for the 
known signal case in the presence of additive background and 
white noise components. In IR detection, the target ampli- 
tude and position are usually not known completely. The 
threshold detector is an extension of Harger' s work which 
accounts for the unknown amplitude and position. 

The form of the optimal detector under the Neyman-Pearson 
criterion will be derived for multiple observations. It is 
assumed that the image plane and object plane remain parallel and 
stationary with respect to each other during the observations. 

The detection of the target under the Neyman-Pearson 
criterion becomes a problem of hypothesis testing. The 
task of the detector is one of choosing between two 
hypotheses, H~ that only noise is present and that in 
addition to the noise there is a target present. The design 
of the detector is one that permits the correct choice of 
hypothesis (a detection) with maximum probability of 
detection, Qd , in a fixed probability, Qp^, of choosing 



20 



H-^ when Hq is true (a false alarm). The structure of the 
Neyman-Pearson criterion requires forming a likelihood ratio, 
Ajq, which is compared to a known threshold, V , specified 
for a given false alarm probability. If the likelihood 
ratio exceeds the threshold, hypothesis is assumed true 
and if not, hypothesis Hq is assumed true. The logic process 
may be written 

«1 

A 10 < V (3.1) 

H o 

where the threshold, V , is derived from the expression 

oo 

Q fa = / P (Z I H 1 ) dZ (3.2) 

V t 

and Z is a sufficient statistic of the received data. 

A. DETECTION OF A KNOWN SIGNAL 

The data for processing is a set of M functions of time 
(Z n (t), n=l,...,M} where each function is the output of the 
preamp during one nutation (0,T). T is the nutation period. 
The output may be represented as 

z n (t) = S n (t) + B(t) +N n (t); 0 < t<T, n=l , . . . ,M (3.3) 

S (t) represents the output due to a target. The component, 
B(t), represents an additive background and is a sample 
function of a random field (B) . The relative size of the 
target in the field of view is assumed sufficiently small 
that the background is additive while B(t) is assumed to be 
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the sample function during each of the M nutations observed. 
N (t) represents the additive noise component due to the 
internal noise of the detector-preamp combination. It is a 



sample function of a white Gaussian random process (N ) of 

2 

zero mean and power spectral density N/2(V /H^) . 

The two hypotheses for the likelihood ratio may be 
written as 

Hr <Z n (t ) = S n (t ) + B(t ) + N n (t ) ; n=l , . . . ,M} 

(3.4) 

H Q : {Z n (t ) = B(t ) + N n (t ) ; n=l,...,M 

Denote the likelihood ratio as A^Q({Z n }) where the like- 
lihood ratio is defined as 

P({Z }|H,) 

A 10^ Z n 1 '* " p ft Z n > | H 0 ) ( 3 -5) 



The likelihood ratio is complicated by the fact Z n and 
Z m are not independent because of a common noise component, 
B(t) . To simplify the derivation Harger introduced an 
auxiliary hypothesis 

H 2 • (Z n = N n ; n= 1 , . . . ,M} . (3.6) 



Then the likelihood ratio, A^g> can be computed using 
the chain rule for likelihood ratios: 



A 
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p({Z n )|H 1 ) p((Z n }|H 1 )/p({Z n }|H 2 ) 

P(^ n >|H 0 ) = PC{Z n >|H 0 )/p(lZ n T|H^ 



A 12' /A 02 
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If the sample functions for the background, B, were non- 
random, Z^(t) and (t) , i f j, would be mutually indepen- 
dent because of white noise. Accordingly, if B were fixed, 

+- 

one can write the conditional likelihood for the n w 
observation ({Z }|B). Moreover, the likelihood ratio 

for M observations would be the product of the likelihood 
ratios for each single observation, 

M 

A 12 C{Z n } l?3 ■ ", <2 

n= 1 



The conditioning is then removed by averaging over the 
ensemble of backgrounds 

" E B[ A 12( { ?n>l« ] ' 

Since the noise is white Gaussian, the conditional 

likelihood ratio is that for detecting a known signal 

S +B in white Gaussian noise and has the well known form 
n n 

( [Ref. 15] , page 253) 



A 



n 

12 ’ 






exp 





Z n (t) (S n (t) + B (t) ) dt 



1 M T 2 

- - 2 / (S (t) + B (t) J dt 

N n=l 0 n 



which may be rewritten as 



A 



n 

12 



ay |b) 



exp 



, M T 

£ E / Z (t) S (t) dt 
N , n n K J n v 
n=l 0 



1 M T 2 

77 B / S‘(t) dt 

N n= 1 0 n 



exp a (b) 



(3.14) 
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where 



a( b ) - Z / (Z„(t) 

n=l 0 

To find A 10 ({Z n }), 
where 



+ S n (t)) B (t) dt - 



one must evaluate 



T 

~ / B 2 (t)dt. (3.15) 
0 

E B[ A \ 2 ( {Z n } l B) J 



E Bt An i2 ({Z n } l B ^ 



? M T 

exp (jqr z / Z (t) S (t) dt 
n=l 0 

(3.16) 

1 M T 2 

' N Z A S n * e b ^ exp a ( b ^’ 

n=l 0 



Assuming the background is Gaussian the expectation can be 
evaluated by representing the background in a Karhunen- 
Loeve (K-L) expansion 

B(t) = Z b k b k (t) » 0 £ b 1 T ( M7) 

k 



where <f>j,(t) is the eigenfunction corresponding to eigenva] le 
of the homogeneous integral equation 

T 

k k <f> k (t) = A Rg(t,u) <J> k (u)du , 0 £ t < T. (3.18) 

Rg(t,u) is the covariance function of background 

Rg (t ,u) = E{[B(t) - m fi (t)][B(u) - m B (u)]} (3.19) 



where nig(t) is the mean of the background. 
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Then Eg(expa(b)) can be calculated in a straightforward 
manner. After the lengthy calculation shown in Appendix C 

In Eg {exp a(b)} = - j Z In (1 + ^ A^) 

k 0 



- T M T M 

* p it £ (Z n (t) - S (t)) / du E (Z.M-S 1 (»)),(t,») 
0 n=l " 11 0 m-1 m m 



T M 

i / dt ( nig ( t ) - l Z (Z (t)-S n (t))} h(t) 
0 n=l 



(3.20) 



where q(t,u) satisfies the integral equation 



_N 

2M 



q(t,u) + / Rg(t,x)q(r,u)dT = R R (t,u) 0 u, t < T (3.21) 



0 



B 



and h(t) satisfies 



_N 

2M 



h(t) + / Rg(t,t)h( T )dT = m R (t) 0 < t £ T 



0 



B 



(3.22) 



q(t,u) is the impulse response of the filter that gives 

A 

the minimum mean square error estimate B of the background 
when the mean is zero [Ref. 15]. The additive white noise 
is reduced in amplitude by 1/M. Using the K-L representa- 
tion, one can also write 



q(t,u) 



Z 

k 



_N 

2M 

N 

2M 



l k 



4> k (t) 4>{(u) 



0 < u , t < T 



(3.23) 



Using equations (2.9), (3.10), (3.14) and (3.20) and 



25 



observing that 



A 02 ({Z n }) “ fl 12« Z n 1 >l S n ' 0 C 3 - 24 > 

it can be shown (see Appendix C) that the likelihood function 
for hypotheses and is 

M T 

A in ({Z }) = exp { Z / Z (t) q (t) dt + KM (3.25) 

iU n n=l 0 11 n 

where 

q n (t) = ^ [S n (t) ' f z S £ (u)q(t,u)du] (3.25a) 

0 k 1 

and 

„ , = 1 M T 7 , TT M M 

N Z / S n (t)dt + ~ ff Z Z S (t)S (u)q(t,u)dtdu 

n=l 0 n 1NU 00n=l 1=1 n 30 

M T 

- ~ Z / S (t)h(t)dt I>.25b) 

m n=l 0 

The logarithm is a monotonic function. Thus it offer 
no loss of information but greatly simplifies the calcula- 
tions. In the remaining work 

A ’l 0 ({Z n » - in A io (( z n n . (3.26) 

Rewriting equation (5.25), 

A 'lO ({Z n }) = L ({Z n }) + K ’ (3 - 27) 

where 
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M T 

L ( Z ) = Z / Z (t)q (t)dt . 
n n=l 0 n n 



K' does not depend on the input data so it may be included 
in the threshold. The hypotheses test now becomes 



L ( Z 



n 




> 

< 



H, 



Threshold . 



(3.28) 



It is also possible, as Barger has shown, to derive 
tests similar to (3.28) without making the Gaussian back- 
ground assumption. In particular, if the background fluctu- 
ates many times over (0, T) such that a(b) of (3.15) is 
approximately Gaussian, then the test is 



M T x 

Z / Z (t) g (t)dt „ Threshold (3.29) 

, _ n n < 



where 

- Tf ' Rb'Cu) [i ^ S l( u)]du] . 

Since by Mercer's theorem 

oo 

Rg(t,u) = Z * k <f> k (t) 4>J(u), 0 £ t, u £ T (3.30) 

k 

on comparing , g^ and (3.21), it is seen that the above 
model approaches the Gaussian model if 
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N « i 
2M X> m £ x X k ' 



The above condition, which implies a weak background, is 
sufficient but by no means necessary. The two models 
approximate each other whenever the second term of and 
are nearly identical. The point is that while the 
Gaussian background assumption may not fit the reality, 
a processor based on the assumption is probably not far 
from being the optimum. Also, the Gaussian assumption allows 
one to compute the processor performance. 

B. DETECTION OF SIGNAL WITH UNKNOWN PARAMETERS 

The likelihood ratio derived in section A is applicable 
only if the signal is known completely. Such is seldom the 
case in IR target detection and in general the signal is a 
function of some unknown parameters. Typically, the ampli- 
tude, a, of the target and its position r Q in the image plan 
are the unknown parameters, and the form of the signal can 
be rewritten to include these parameters: 



S = S (t;a, r) 
n n^ ’ ’ o' 

Furthermore, by assuming that only point targets are of 
interest and including the assumption of fixed image and 
object planes over the number of nutations of interest, 
the shape of the waveform will be the same during each 
nutation. However, the amplitude is allowed to vary between 
successive nutations. 
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The signal waveform may now be written as 



S ( t • a , r ) = 
n v ’ ’ o J 



a f (t , r ) 
n v ’ o' 



(3.31) 



where a n is the signal amplitude and f(t,r ) is the signal 
shape for a target at r . 

Rewriting hypothesis under these assumptions 
gives 

Hj’: (Z n (t) = a n f(t,r o ) + B(t) + N n (t); n=l,...,M 

0 < t < T) (3.32) 



The likelihood ratio A^q becomes the conditional 
likelihood ratio A^q ({ Z^} | a , r ) obtained by replacing 
S (t) by a f(t,r ) everywhere in (3.25). 

Assuming a priori statistical knowledge of a and r , 
and their independence, the likelihood ratio may be 
averaged over the respective density functions to produce 
an unconditional likelihood ratio 

A 10^ Z n^ = // A 10^ Z n^~ , ^o)P(~)' > (^o) d ^ d ^o ’ (3.33) 

Instead of averaging as above, one could estimate a and r Q 
by the maximum likelihood procedure and substitute the 
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estimated values a and in the conditional likelihood 
ratio and try to maximize the generalized likelihood ratio 

A 

A -y 

AlO ({Z n } I a , r o) • However, this procedure is not followed 



here because the maximum likelihood estimates are difficult 
to ob.tain. 

Since the main interest is in detection of the \\ T eak 
target, the small signa] case is of primary interest. 
Expanding } ' a , r Q ) in powers of signal amplitudes 



M 



A i n I a ,r n ) = 1 + Z a 

n=l 



A in ({Z }|a,rj 



1 0 v 1 "n J 1 ~ ’ * 0 J “n 3a “10 v liJ n J 1 “ * 1 Cr 1 a=0 

n 



(3.34) 



1 M M a 2 

2 a n a k 3a 3a, A 10 ({Z n } I ~ ,r 0^ a = 0 

n=l k=l n k 



Integrating (2.34) as in (3.33) term by term and 
retaining only the first non- zero term that depends on the 
data provides a statistic which is optimum for weak signals 
and is called a threshold detector [Ref. 4]. 

The evaluation is made easier by the fact 



3 

3 a. 
k 




uVh-V 



a=0 






which results from 

A({Z n )|0,r 0 ) = 1 . 



(3.35) 
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The first non-zero term becomes 




- ^ / z n (t)dt / q(t,u)f (u,r Q )du 



T T 



0 " 0 



T 



( 3 . 36 ) 



/ f (t,r n )h(t)dt 
0 u 



Performing the integration over r^, 



M 



T 



T 



A 10^ Z n^ S 1 + N 2 *n' f Z n ( ' t ' ) “ / f ( u ) du ] 



where the bar denotes the averaged quantities. The first 
and last terms are independent of the data and may be 
included in the threshold value. Moreover, if the ampli- 
tudes are identically distributed, as is generally the case, 
the threshold detector performs the test 



n=l " 0 



T 



-/ f(t) h(t) dt 
0 



( 3 . 37 ) 



M T 

W = Z / Z (t) q (t ) dt 
n-1 0 n 




( 3 . 38 ) 



where W is the threshold and 



q(t) = ^ [f(t) - / q ( t , u ) f(u) du] . 



0 



( 3 . 39 ) 
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Thus the threshold detector can be mechanized either as a 
correlator or matched filter followed by an accumulator. 

The filter function q(t) may be represented in a series 
expansion in terms of the eigenfunctions of Rg(t,u) and the 
coefficients of the signal expansion. Expanding f(t) in 
terms of the eigenfunctions: 



00 

f(t) = z f k 4> k (t) . 
k 

Substituting (3.40) in (3.39) and using (3.23) 



q(t) 



_N 

2M 



N 

2M 



+ 



*k (t) 



(3.40) 



(3.41) 



C. EVALUATION OF PERFORMANCE 

The Gaussian assumption was made in deriving the 
threshold statistic; therefore, the threshold statistic 
itself is Gaussian. This means that the probability of 
detection, Qd , and probability of false alarm, Q^, are 
uniquely determined by second-order statistics, e.g., 
the means and variances of W under hypotheses Hq and H^. 
Recall that the threshold statistic is 



M T 

W = Z / Z (t) q(t) dt. 
n=l 0 n 



(3.42) 



The means and variances of W are: 
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M T 

E(W|H ) = E / q(t) E (B (t) + N (t)) dt 
n-1 0 



M T_ 

= 2 / q(t) m (1.) dt 

n=l 0 B 



M T 

E(W|H ) = E / q(t) E [ a f(t,r 0 ) + B(t) + N (t) ] dt 
n=l 0 



M T 

= E(W|H ) + E / f(t) q (t) dt 
u n=l 0 



VAR (W|H Q ) = VAR (W 1 H 1 D = VAR 



= / q^(t) dt + / q ( t) / q(u) R R (t,u) dudt 

^ 0 0 0 



(3.43) 



Letting p (»1FL) be the Gaussian density function of 
W under FL which is specified by the equations above, then 

oo 

Q fa = / p w (xlH 0 )dx (3.44) 

K t 

/ 

and 



Q H = / P w (*lH,)dx (3.45) 

d Wt w 1 



where W is the threshold voltage. 

Equations (3.43) and (3.44) can be solved to yield 
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Q fa " M x > 



(3.46) 



and 



Q d = * c (x-d) 



(3.47) 



where 

W t - E[W|H q ] 
/VAR 



d is the equivalent signal- to-noise ratio 

E(W|H ] 

d = — 

/VAR 



M T. 

Z a / £(t) q(t) dt 
n=l n 0 



(3.48) 



/ ^ p ^ 

Jyr f q 2 (t) dt + M 2 / q (t) / q(u) R R (t,u) dudt 
V L 0 0 0 

and 4> c (g) is the complementary error function of the kind 



1 -v 2 / 2 

4>,(g) = / e dv 

d /H g 



(3.49) 



By setting a desired probability of false alarm, x 
may be determined from equation (3.46). The probability 
of detection is then determined by solving for d and using 
equation (3.47). 

By substituting (3.41) and (5.42), and using (3.23), 
the numerator of d is 
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2 



M T 

Z a / f(t)q(t)dt = 
n= 1 n 0 




(3.50a) 



while the denominator of d (see Appendix D) is 




(3.50b) 



Therefore the effective signal-to-noise ratio is 




(3.51) 



Knowledge of d along with (3.46) and (3.47) completely 
specifies the performance of the threshold detector. 

The equivalent signal-to-noise ratio for a target whose 
amplitude and position are known is easily shown to be 



By letting the background be zero it will now be shown 
that d' reduces to the equivalent signal-to-noise ratio for 
a known signal in white Gaussian noise which is well known 
[Ref. 4]. 

The energy dissipated during one observation in a resis- 
tor of 1 ohm if the signal a n f(t,r Q ) is the voltage across 
that resistor is 




2 



(3.52) 
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T oo a> 

Sn = a „ ' 3 f k(V*k(‘5 5 £ J Crjqwit 

O K 36 

oo 

• a n £ l f k( ? o>l • C 3 -53) 

The voltage resulting from M observations then is 



/E 





1 M 
M n=l 



n 



E 

k 



Mv" 




(3.54)' 



Substituting (3.54) in (3.52) with A ^ = 0, 



d' = /E 





/M 






which shows the familiar result that detectability depends 
only upon the signal energy and the signal-to-noise ratio 
increases as the square root of the number of observations. 
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IV. THRESHOLD DETECTOR FOR NOTATING SYSTEM 



In an infrared detection system the output voltage to 
be processed is obtained from a detector-preamplifier com- 
bination. The irradiance on the detector in W/cm is con- 
verted to volts with a linear scale factor. The voltage 
output can be written as 

v(t,r Q ) = KH(t,r 0 ) (4.1) 

where H(t,rQ) is the irradiance on the detector from a target 
at position Tq , and K is a scale factor in (V/W/cm 2 ). Also 
the covariance function may be written as 

R vB (t) = E [ KB ( t + x) KR (t) ] = K 2 R b (t), (4.2) 

where B(t) is the irradiance on the detector due to a partic- 
ular background scene. The output voltage due to this see: i 
is 

v B = KB (t) (4. ) 

The problem now is to relate the detector-preamplifier 
output to the equations derived for the threshold detector. 

It was shown in Part III that an optimum processor for 
weak signals was a threshold detector that performed a 
certain linear operation on the data. The linear filter was 
specified in terms of the eigenvalues and eigenfunctions of 
the integral equation (3.18), where the kernal R B (t,u) is 
the covariance function for the background. 
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In Part II, the covariance function was noted to be 
periodic. Rewriting (4.2) in terms of the output voltage 



7 jkw (t-u) 

R vR (t,u) = k 2 e 

k 



(4.4) 



Substituting (4.4) into integral equation (3.18), the solu- 
tion for the eigenvalues and eigenfunctions is trivial with 



and 






jku. t 
e 

/~T~ 



x 



k 



K 2 B k T 



(4.5) 



where T is the nutation period. 

The signal coefficients were also noted to be periodic 
and could be expanded in a Fourier series 

00 jkw t 

v ( t ) = K E H r (r ) e 0 . (4.6) 

k 



Recognizing the similarity between equations (4.6) and (3.41) 
the signal coefficients for a known signal which is expanded 
using the basis set above become 

VV ' f 4 - 7 ’ 

For the threshold detector, the signal is averaged with 
respect to an a priori density function for r R . Averaging 
can be done with the coefficients to give 
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(4.8) 



f, = / p (r ) f . (? ) dr 
k 1 v o' k v o' o 

= K/T / I 1 , (r )p(r )dr 
k v o' 1 o o 



■ K/T H k ■ 



Using equation (2.10), the coefficients b' k can be written 
as 

d 2 k (4.9) 

The transform of the radiance distribution from a target, 
N'(k,r ) is a function of the target radiance, R(k), the 
position of the target, r , and the point spread function of 
the optics, F (k) . Rewriting (4.9), 



H, 



= j n Jdr Q p(r o ) /N ' (k ,r Q ) t* (k) e"-* 11 ^ 



n 



2ttp -Jk 2 + k 2 



y 



H k = j n jd 2 kR(k) t * (k) J R ( 2 7T p yj k 2 +k 2 e~> n * 



r - i 2 tt k , r 

Jdr 0 p(rje 



(4.10) 



If the density function, p(r Q ), is assumed Gaussian with a 
standard deviation p Q , then the averaged coefficients become 



Hr = j 



n / r oo' 



-2TT P 2 (k 2 + k 2 ) 
l o v x v' 



(k)e 



-jn<(» 



n 



2TTD^k 2 + k 2 d 2 k. (4.11) 



The Gaussian assumption on r Q is reasonable because r Q 
often represents a pointing error from some designated tar- 
get position. Then p Q describes the accuracy of the pointing 
mechanism. Substituting (4.8) into (3.42), the filter func- 
tion for the threshold detector becomes 






q(t) = K Z H k v / T 



k 



1 



\ 



i * f K 2 e k T ; 



i k V 



/T 



(4.12) 



39 



Rewriting the equivalent signal -to-noise ratio in terms 



of equations (4.8) and (4.12) one has 



d = 



M 

Z a /T/M 
M i n 

n=l 

/n7T‘ 



Z [ H, 
k 1 



1 + 



MTK 2 3k 



n7T 



(4.13) 



It will be shown in Part V that 8, may be represented as 



e v = c R V Q k 



2 . 



where cr R is the amplitude of the background noise correla- 

2 

tion function, fi is the instantaneous field-of-view of the 
detector and is an integral function of the optical and 
detector parameters and correlation length of the background 
radiance . 



Define the following terms 



Signal - to -noise ratio (SNR) = 



1 M . _ 
K fr Z a /T 
M , n 
n=l 



✓N/2 



(4.14a) 



Background - to -noi se ratio (BNR) = 



OgQK/T 



/N/2 



Using (4.14a) and (4.14b), d may be rewritten as 

2 1 



d = SNR /M 



£ |H k | 9 

k 1 + MQ, BNR“ 



(4/ 4b) 



(4.15) 



Noise equivalent irradiance (NET) may be defined at the 
detector-preamplifier output to be 
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K-NEI 



(4.16) 



■ 7 ? - 2b w - 

where B is the detector-preamplifier bandwidth, 
w 

The SNR and BNR may now be written in terms of the NEI as 

M SNR = — = SNR' /NO- 
NE I 

where NTB is number - time -bandwidth product and 

Q D 

/M BNR = rnfr /NTB = BNR' /NTB. 

NEI 

Expressing d in terms of NEI one has 



d = SNR’ /NTB E |H 



k' 



1+0, NTB BNR 
K 



,2 



(4.17) 
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V. SOLUTION FOR RECTANGULAR DETECTOR 



The equations in Part II for the detector -preamp 1 i f ier 
output voltage and the correlation function of the background 
must be solved to obtain the filter function and specify the 
detection system performance. 

The two general equations obtained in Part II are 




t* h 

where H (r ) is the n r Fourier series coefficient for the 
n ^ o' 

radiation on- the spatial filter and 



3 



n 



/|tooi 2 



F (k) 
o ~ ' 




where B is the n 1 "^ 1 coefficient of the background correlation 
function at the spatial filter. 

A rectangular shaped detector was chosen for the calcu- 
lations because it exhibited the simplest form of the 
equations. For circular nutation the filter is represented 
graphically in Figure 5.1. 
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Figure 5.1. Nutation with Rectangular Detector. 



If the detector function is given by 



t (x 1 ,y ' ) = 



.. x -w/2 < x' < x +w/2 
1 o'— — o 

y 0 -h/2 < y ’ < y 0 + ^/2 

0 otherwise 



then its transform is 



T( W " 



- j 2 tt fk x +k v ) 

J v x o y o J 

- = sin (irk w) sin (irk h) 

"Vy Y 



so that 



r (k) | 



sin (nk^w) 

272 
u k x 



sin (irk y h) 



(5.1) 



The point spread function of the optics is assumed to 
be Gaussian: 
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f (r) = 
o ~ J 



2a' 



2ira ' 



Its transform is 



F (k) = e 

o ~ J 



, 22 ,, 2 , 2 . 
-2tt a (k x +k y ) 



(5.2) 



The correlation of the background radiance is assumed to 
be of the form 



4, B (x,y) = o 2 e "°l x l e "8 |y | 



(5 . 3a) 



where a * and 3 * are the correlation lengths in the x- and 
y- directions, respectively. The Wiener spectrum is the 
Fourier transform of 4*^: 

23 



2 2a 

*B« * °B ~ 



a‘ * (27rk x ) 2 8 2 ♦ (2irk y ) 2 



(5.3b) 



Substituting (5.1), (5.2) and (5.3b) into (2.10), the 
radiance from a point target becomes 






■ ■"/■ 



j 2 tt (k r +k r ) -2ir 2 a 2 (k 2 +k 2 ) 

e J v x x y y J e 1 x y J 



i 2 tt (k x +k y ) simik w sinmk h 
J K x o y'o' x y 

Fk tt k 



y 



-jn tan " 1 (k /k ) 



y x J ( 2 tt p /k 2 + k 2 )dk dk . (5.4) 

nV vxyixy 



Letting 



u x = y/r 27rok x 



u y = /2iroky 
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equation (5.4) becomes 



_ -nr -(u 2 + u 2 ) 

i- c x y 



H (r ) = 2— / e 
n^ o' 21 



-j — f~(r -x )u + (r -y )u "j 
. n ° L x o' x ^ y 'o' yj 



; jn /rot 



n\ o x/V u y du x du y ( 5 ' 5 > 



If r Q is random, the radiance from a target may be 
averaged with respect to the density function p(r Q ) as shown 
in Part IV. In particular, if r Q is Gaussian distributed, 
with standard deviation p Q the averaged coefficients become 



x . n f 

n J 



-2h*y* (k 2 +k 2 ) j 2 n (k x +k y ) 



x y 



x o y' o' 



-1 



sinuk w sinuk h -in tan (k /k ) 
x y J K y x' 



where 



'n ( 2 " p /?? 



k 2 +k 2 I dk dk 
x y x y 



p + a 
K o 



Letting 



u ' = /Tiryk 
x ' x 



u ' = /2"iTYk 

y y 



( 5 . ( i 



equation (5.6) becomes 
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'2 '2, , /2' , 



li = CJe' (U ^' U ^ e’ jn (u "x /u y ) 



n 2 

TT 



n 



— f - /u' 2 + J 2 I du ' du ' 






Y 



x y I x y 



(5.7) 



And finally substituting equations (5.1), (5.2) and 
(5.3b) into equation (2.22), the correlation coefficients 
become 



n 






- 4tt 2 o 2 (k 2 + k 2 ) sin 2 (-rrk^w) sin 2 (irk^h) 



(Kk x ) 



O k y) 



2a 



2 B 



B 2 ± m i ^2 .2 , , ,2 

a + (2irk x ) 6 + (2irk y ) 






Letting 



dk dk 
x v 



(5.8) 



u" = 2Tipk 
x K x 



Uy = 2rpk y 



equation (5.8) becomes 



o 2 f2 2 a6a 2 r -u“-u 
6 = e x y 

n u 2 J 



w 



2 - h 



2 "2 /Sin — — u" \ /sm u 

■■ X \/ /— _ V 



/2a 



w 



u’ 



o' X 



/2a 



/2a 



U i 

' /2o y 1 



1 



(aa) 2 +u 2 (6a) _ y 



1 t 2 ( P /”2 , "2 \ 

i 2 + u 2 n u y j 



u^+u" idu"du" (5.9) 

x y I x y ^ ' 



where 



n 2 2 , 2 
ft = w h . 



(5.9) 
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The integrals in (5.5), (5.7) and (5.9) were difficult 
to evaluate. Although the components of the integrand are 
relatively simple functions, the argument of the Bessel 
function prohibited separating the integral into two one- 
dimensional integrals. No closed form solutions were found 
so the integrals were evaluated numerically. The numerical 
method used is based on w r ork by Pierce [9] who applied 
Gaussian quadrature formulas to two dimensional integration 
by integrating over a planar annulus in the (x,y) plane. 

Gaussian quadrature formulas are a means of evaluating 
the integral by summing weighted values of the integrand at 
specific points. For the one dimensional case, 

b b 

f g(x)dx = / w(x)f(x)dx (5.10) 

a a 

where g(x) is the integrand to be integrated and w r (x) is a 
weight function for which the specific integration was 
derived. For well behaved functions the integral may be 
evaluated as 

b m 

/ g(x)dx = E A.f(x.) + error (5.11) 

a i=l 1 1 

where the x^'s are the m zeros of the m^ polynomial 

m 

P (x) = II (x-x.) 
m i = l 1 

of the set of polynomials mutually orthogonal on the inter- 
val [a,b] with respect to the weight function w(x) . The 
weights are given hy 
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(5.13) 



b 

= / w (x) (x) dx 

a 

where (x) = P m (x) / [ (x -x^ (x^ ] is the Lagrange inter- 

polation coefficient. 

For a weight function w(x) and interval [a,b] the poly- 
nomial P (x) and its zeros x. and weight factors A. need be 
computed only once. For a number of weight functions and 
intervals the set of orthogonal polynomials is known. Stroud 
and Secrest [13] give x^'s and A^'s for a variety of weight 
functions and internals. The degree of the formula deter- 
mines the number of x^'s and A^ ' s . A 2M-1 degree formula 
will have M points and is exact for polynomial integrands of 
degree 2M-1 or less. 

Pierce applied Gaussian quadrature integration to two 
dimensions where the solution is of the form 



The integration is over an annulus in the x,y plane wi h 
inner radius R and outer radius 1. Rewriting the integral 
in polar coordinates 




1 2 TT 

/ / 

R o 



rG (r , 0) d0dr -EE D..G(r.,0.) 

i j J 3 



(5.14) 



Pierce showed the summation could be rewritten as 



4(m+l) m+1 

E E A. B • G (r • cos0-,r. sin0.) . 

i = 1 j = l -’ J J 



(5.15) 
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where and Bj are the weights for the radial and angular 
directions respectively and 4M+3 is the degree of accuracy. 

For an arbitrary annulus the formula may be obtained by 
rewriting equation (5.14) 

2 7T r 2 

I = / / f(r,0) rdrde (5.16) 

0 r, 

and letting 



P 



I 



r 




2tt 

/ 

o 



1 



f 




pf (r 2 P ,6) 



dpd0 . 



The approximate solution can be derived as 



I 



4 (m+ 1) 

Z 

i= 1 



m+ 1 
Z 

j = l 



A.B.f 
i J 





where 



(5.17) 



(5.18) 



A^ = 2tt/ 4 (m+ 1 ) 

and 

B j = w j( r 2 ' r i)/ 2 ’ w j = weight 

M+l is the order of the orthogonal polynomial, in this case 
the Legendre polynomial on the interval (0,1) andcr 's are 
the zeros of this polynomial. This type of formula is known 
as a spherical product Gauss formula [12]. 

The annuli were picked at radii corresponding to the 
zeros of J q (z) because J q (z) is the fastest damping component 
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of the integrands. One additional annulus inside the first 

zero was found necessary to account for the — term when 

a+x Z 

a is small and x (and y) approach zero. The distances 

between annuli outside the first five to eight zeros were 

found not to be critical and thus spaced arbitrarily. 

The Gaussian Quadrature formula used was a Gauss - Legendre 

24-point formula as listed in [13] thus giving a degree of 

accuracy of 99. The integration is limited to only the first 

quadrant as mentioned below. Thus for a 24 -point formula, 

2 

(M+l) or 576 points per annulus were used to evaluate the 
integral . 

The expression for $ n is easily seen to be an even func- 
tion both radially and about both planer axes, therefore, 
only integration over the first quadrant was necessary. 

With a little more difficulty, F (r ) and H can be shown 
(see Appendix E) to exhibit the same property provided one 
expression is used for even n and another for odd n. By 
calculating only over the first quadrant computer time was 
considerably shortened. 

In order to solve for coefficients of arbitrary order, 
the Bessel functions of that order must first be obtained. 

The available library subroutines were found inadequate for 
this use for two reasons. First, these routines have an 
upper limit of 100 on the order of the Bessel function that 
can be computed and second, the routine had to be called for 
each individual order. 
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The recursive relation 



2n 






(5.19) 



provides a basis for generating an array of Bessel functions 
of various orders. However for accuracy J n+ ^(vj and 
must be known. 

The most accurate method to numerically calculate Bessel 
functions of various orders and arguments was found to be a 
uniform asymptotic expansion involving Airy functions [1] : 

’ a^CA) 



J n (nz) - 



4 A 
1 - z 



18 lA l( n 2/3 A) 



n 



T7T 



k=0 n 



TF 



A.'(n 2/3 A) co b k (A) 
573 .. 2k 



n 



k=0 n 



(5.20) 



where A. and Aj are two types of Airy functions and 



I A 3/2 = in 



1 + 



“ z 



- n/i-z 2 



< 1 



a 0 (A) = 1 



b (A) = -5/48 A 2 + — 
0 /A 



24(l-z 2 ) ^ /2 



1 



8C1-Z D 2 



< 1 



a k (A) << a Q (A) 
b k (A) << b Q (A ) 



k > 1 



(5.21) 



The Airy functions are calculated using the equations 



A i (x) = j x"" 2 e~ 5 f (-6) 



and 



A[(x) = - j e' 6 g(-6) 
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where 




(5 .22) 



and 




The functions f(-6) and g(-o) are tabulated and linear 
interpolation is used to determine their values. A is guar- 
anteed to always be positive if n is selected to be larger 
than the argument. 

By picking an order much higher than the argument, v, 
and J using the asymptotic expansion, the recursive 

relationship (5.19) may be used to generate an array of 
successive orders of the Bessel function down to and includ- 
ing J Q . However, error builds up rapidly using this method. 
A second relation 



Multiplying each Bessel function value generated in the 
recursive relation by C resulted in very accurate values 
being obtained. The computer program used to solve the 
equations is listed as Appendix H. 

To evaluate the numerical results a check was used which 
summed the coefficients and compared this summation to an 
analytic expression for a sum that could easily be solved 
numerically . 



1 = J Q (y) + 2J 2 (y) + 2J 4 (y) + 



(5.23) 



may be used to generate a normalizing factor 
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For a particular phase angle of the nutation cycle, w 0 t, 
the radiance function for a target is 



H O 0 t) = I H k (r 0 )e 






(5.25) 



For a point target and a rectangular, circularly nutat- 
ing detector 

2 2 2 

- 2 tt o k siniTwk j 2irk (pcosco t - (r -x )) 
H(w 0 t) = /e x — F — x - e x ° x ° dk. 



fe 



2 2 2 

‘ 2tt a k sin^hk y j 2irk (psina) t - (r -y )) 

■ ' „k.. ' e dk y 



y 



(5.26) 



which has the solution (see Appendix F) 



f 



F l(r x ,x 0 ,w) - V r x-V wl 



- E 2 ( r y ,y 0 ,-h) 



E 2 (r y ,y o h ) 



where 



E 1 (r,z,a) = Erf 



r- z - pcosoj t 
o_ a_ 



/2a 



2 /2c 



E ? (r,z,a) = Erf 



rr-z-psimo t 
o_ a_ 



/2a 



2 /2a _ 



and Erf is the error function 



? x _ 2 

Erf(x) = — / e u du . 

/F 0 



(5.27) 



By evaluating the analytic expression at the phase in- 



stant (w t) that the detector crosses the target and comparing 
v o 
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it against equation (5.27) for an appreciably high n, a 
reasonable check was made on the accuracy of the individual 
signal coefficients. 

For the coefficients averaged over r , the check was 

made at to t = 0 where 
o 



H (o) - Z H, (5.28) 

k k 



and the analytic expression is 

11 Co) = J [ E l (r x- k 0 - w) ■ E l (r x- x o’- w )][ E 3C r y>'o h C 

E 3 (r y> y o' h 3 

where 



E 3 (r,z,a) = E 2 (r,z,a)l t=Q (5.29) 

1 o 

The rms background voltage was used to check the coef- 
ficients in the correlation function. The mean square vol - 
age is 



e 2 = E{B 2 (t)} = Rg (0) = Z S k . (5.50) 

h 

Appendix B shows that the correlation function may be written 
as 



R B (T) = / I T (k) | 2 W’(k) 

Substituting (5.1), (5.2), and 
ing at x = 0 



J 



o 



4iTpk 

i 



% T 
sin - T - 




(16B) 



(5.36) in (16B) and evaluat- 
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O A 2 2, 2 

sin (irk vv) 71 a x 

e B = °B 2a ; r *2' "2 ~ 7 72 dk x 

(Trk x ) a + (27ik x ) 



o a 2 2, 2 

• 2., - 4 7T a k 

sm (nk h) e y 

1 o D 2B / X — T dk 

B n 2 „2 . .. 2 y 



C^y) 



B + (27Tk y ) 



(5.31) 



It can be shown (see Appendix F) that an analytic expres 
sion is 



e 2 = 4aBa 2 R(w;a)R(h;B) 



where 



R(y,a) = 



-u 


1 - Erfc l^r 


2 a 


2 

a 


1 a 


/rfy 



1 - e 



- riL.-) ^ 
l 2a J 



2 2 
a a 



2a‘ 



e ay Erfc 


V 

aa + — 




2a 



+ e 



‘ ay Erfc 



ao 



2 ^ - 2Fr£c(a a) 

and Erfc is the complimentary error function 



(5.32) 



00 2 

Erfc(x) = — / e u du . 
/tF x 



(5.33) 



The computer program used to solve equations (5.6), 
(5.7) and (5.9) and calculate the infinite summations is 
listed as Appendix H. 
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VI. NUMERICAL RESULTS 



Numerical results constitute solving for the signal and 
background coefficients and using these coefficients to 
specify the form of the filter to be used in the threshold 
detector as well as to calculate the probability of detection 
as a function of s ignal - to-noise ratio (SNR), background- to- 
noise ratio (BNR) and number of nutations (M) used to make 
the decision. The probability of detection is determined 
for the averaged equivalent signal -to-noise ratio and for 
specific point targets. 

Many different values for system parameters and target 

locations were used to compute the coefficients. A typical 

set of parameters i s ^ : 

nutation radius = 15 
width of detector = 30 
height of detector = 1 
blur circle standard deviation = .5 
background correlation length x-airection = 20 
background correlation length y-direction = 20 
position of detector in detector coordinate system 
x direction = 15 
y direction = 0 

target coordinates (if desired) 
x direction = 0 
y direction = 0 

pointing error standard deviation = 7.5 
A detector utilizing the parameters above would trace out 
the area in the image plane shown in Figure 6.1. 



All units are in milliradians . 
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Figure 6.1. Scanned Area of Nutating Detector. 

This type scan was chosen because it crossed the center 
of the coordinate system which would be the center of the 
target blur circle for zero pointing error. Also, in the 
case of a random pointing error in the absence of bias ze d 
mean is the center of the coordinates. This is important 
because the threshold detector design is based on a Gaussian 
pointing error with zero mean. 

A sample calculation of 121 coefficients using the para- 
meters above is listed in Appendix G. The equations shown 
in Part V were used to check the accuracy of the coefficients. 
Table 6.1 lists typical values for the summation of 121 coef- 
ficients for various sets of parameters and the calculated 
summation values. 
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TABLE 6.1 



Calculated Sum 



Actual Sum 



Signal (point target) 



.6826894 

.6826895 

.5438002 



.6826993 
.6826994 
. 5438087 



Signal (averaged) 



. 2990664D-6 
. 2 32 81 38D- 4 
. 2651595D-1 



. 2991232D-6 
. 2328142D-4 
. 2651595P-1 



Background 



.6013477 

.2347420 

.8335525 



.6012230 

.2347420 

.8319764 



Besides the general shapes of the background and signal 
spectra, variations in the spectra due to detector size, 
background correlation lengths and coordinates of the detec- 
tor are of interest. Curves in the following figures show 
the effects of these variations. 

Figures 6.2 through 6.5 show the envelope of the one 
side of the two-sided background power spectral density 



where w is normalized to 2 tt . The curves of Figures 6.2 and 



6.3 illustrate the effect of varying the background correla- 
tion lengths. In Figure 6.2 the correlation lengths are the 
same in both x- and y-directions . Practical measurements, 
however, indicate that correlation lengths in the horizontal 
(x-) direction may be longer than in the vertical (y-) direc- 
tion. The curves in Figure 6.3 are for this case. 



oo 



S B (u) = T. B k 6(w-ku o ) 



k 



o 
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The curves of Figure 6.4 and 6.5 relate the spectrum to 
changes in size of the detector. In Figure 6.4, only the 
detector with BNR equal to 9 and number of nutations equal 
to 100 are shown in Figures 6.11 through 6.13. Each figure 
shows the envelope of one side of a double sided spectrum. 

Two cases are considered. First, threshold detection systems 
designed for different Gaussian pointing errors but with the 
same size detector are shown in Figures 6.11 and 6.12 and 
second, detection systems designed for the same printing 
error but utilizing different size detectors are shown in 
Figures 6.12 and 6.13. 

The curves shown in Figures 6.14 through 6.16 show the 
matched filter to a particular point target located within 
the scan. Notice each of the filters exhibits a band pass 
characteristic. This is because of the additional noise 
component (background) . The oscillations observed in the 
point target spectra of Figures 6.6 to 6.8 have been sup- 
pressed to provide a look at realizable filters. 

The probability of detection for the threshold detector 
was shown in Part IV to be 



Q , = SNR M 



l 

k 



+ MQ k BNR‘ 



Figures 6.17 through 6.20 show as a function of SNR, 

BNR and M . 

The probabilities of detection for point targets located 
at specific points in the scan have also been plotted. The 
probability of detection for a specific point target is shown as. the 



So 



width of the detector is varied. As width is increased, so 
is the area scanned by the detector. Figure 6.5 shows the 
difference in spectra between a rectangular detector and a 
square detector which have approximately the same surface 
area. It should be obvious however that the rectangular 
detector scans much more area per nutation than a square 
detector of identical surface area. 

The envelopes of signal coefficients (magnitude only) 
for a point target are plotted in Figures 6.6 through 6.8. 
Phase information depends only on the location of the target 
relative to the initial point for the nutation cycle which 
is along the x axis. This has not been plotted. Changes in 
the magnitude due to changes in detector size are shown. 

Figures 6.9 and 6.10 show the envelopes of coefficients 
for the signal that is averaged with respect to a Gaussian 
pointing error. The bandwidth of the averaged signal is 
much less than a point target as expected. The figures she 
the spectra for two different standard deviations of point- 
ing error. 

Matched filters (magnitude only) for the threshold 



= SNR'/M - 





k 1+MQ^BNR 



where 



SNR’ = i 



M 




a 



n 
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Figures 6/21 through 6/24 show as a function of SNR', 

BNR and M. 

The parameters of the nutating optical system have been 
varied to show their effect on 0^ and Q^. 
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Figure 6.2. Background Power Spectral Density vs. Frequency (Hz). 

Detector parameters: RHOIO, SIG=.5, w=10, h=l 
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Figure 6.3. Background Power Spectral Density vs. Frequency (Hz). 

Detector parameters: RHO-15, SIG=.5, w=30, h=l 
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Figure 6.4. Background Power Spectral Density vs. Frequency (Hz). 

Common detector parameters: RH0=15, SIG=.5, h=l 

Background parameters: a=$=.l 
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Figure 6.5. Background Power Spectral Density vs. Frequency (Hz) 

Common detector parameters: RH0=15, SIG=.5 

Background parameters: a=B=.l 
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Figure 6.6. Point Target Spectrum Envelope vs. Frequency (Hz) 

De^--i-or navameters: RH0=15, SIG=.5, w=30, 

h=l, x q = 1 5 , y Q =0 

Target location: R =0, R =0 
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Figure 6.7. Point Target Spectrum Envelope vs. Frequency (Hz). 

Detector parameters: RH0=15, SIG=.5, w=10,h=l, 

x o =15 > >V° 

Target location: R =0, R =0 
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Figure 6.8. Point Target Spectrum Envelope vs. Frequency (Hz) 

Detector parameters: RH0=1S, SIG=.5, w=5.5, 

h=5.5, x Q =15, y o =0 

Target parameters: R =0.0, R =0.0 
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Figure 6.9. Averaged Signal Spectrum vs. Frequency (Hz). 

Detector parameters: RH0=15, SIG=.5, w=30, h 

x =15, y =0 
o ’ ' o 

Std. Dev. of Gaussian Pointing error = 7.5 
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figure 6.12. Filter Spectrum of Threshold Detector vs. Frequency (Hz 

Detector parameters: RH0=1S, SIG=.5, w=30, h=l, 

x =15 , y =0 
o ’ 1 o 

Std. Dev. of Gaussian pointing error = 7.5 
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Figure 6.13. Filter Spectrum for Threshold Detector vs. Frequency (Hz) 
Detector parameters: RH0=1S, SIG=.5, w=5.5, h=5.5, 

x o =15 > >V° 

Std . Dev. of Gaussian pointing error - 5.0 
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Figure 6.14. Fil’ter Envelope (Matched to Point Target) vs. 
Frequency (Hz) . 

Dete^c. '.uiameters: RH0=15, SIG=.S, w=30, h 

x o =15 > y o =0 

Target location: R =0, R v =0 
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Figure 6.19. Probability of Detection vs. Signal -to-Noise Ratio. 

Detector Parameters: RH0=15, SIG=.5, w=5 5 h=5 5 

V= 15 ’ >V° ' ’ 

Std . Dev. of Gaussian pointing error = 5.0 



o 






• 






o 




CD 


•H 

4 




a 


e3 




O 






_ * — 1 


<D 


#N 




CO 


1 ( 




•H 


II 




O 


X 


CD 


i 




« 


O 


o 


_cn 


4 


hO 




l 


II 




T3 


£ 




c 






3 


#v 


CD 


o 


LO 


% 


4 


• 


-CO 


bo 

rX 


II 

O o 




u 


l-H II 




03 


CO o 




PQ 


X 


o 


. 


#v 

LO 


' • 


CO 


r— 1 LO 


;_in- 


> 


11 rH 

O II 




P 


X o 




o 


x 








! ° 


4 






O 


• • 


Lai 


<D 

4 


to 

4 




<L> 


<U 




Q 


4 






<D 


o 


4-1 


& 


O 


o3 


* 




4 


_LO 




o3 




4 


Ph 




•H 






rH 


4 




•H 


O 


CD 




*4 


03 


U 


» 




<U 


_ vi 1 


O 


4 




4. 


<U 




C4 


Q 


CD 


# 






O 




-_.ro 


CVJ 






\o 






<u 




o 


4. 

3 




a 


bO 




_C\J 


•H 






Ph 





LO 

II 



cn 


CO 


rv. 


CD 


UP 




OP 




r-| 


o 


O 


o 


CD 


o 


ci 


o 


o 


CD 


O 


o 



80 



Std. Dev. of Gaussian pointing error 
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Figure 6.22. Probability of Detection vs. Signal- to-Noise Ratio. 

Detector parameters: RH0=15, SIG=.5, w=10,h=l, 

x ° =15 > y 0 =o 

■ Std. Dev. of Gaussian pointing error = 5 
Target location: R =0, R =0 



SNR = 100 



D 
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Figure 6/23. Probability of detection vs. Background- 
to-Noise Ratio. 

Detector parameters: RH0=15, SIG=.5, w=30 

h=l,x o =15, y Q =0 

Std. Dev. of Gaussian pointing error = 7.5 
Target location = R =0, R =0 
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Std. Dev. of Gaussian pointing error 
Target location = R =0, R =0 



VII. CONCLUSIONS 



The problem investigated in this thesis was to determine 
an optimum processor to be used in detecting targets with an 
infrared nutating system. The nature of the system divided 
the problem into two areas: optimizing a spatial processor 

according to size and shape and optimizing a temporal pro- 
cessor that takes the output time varying voltage and decides 
if a target is actually present. 

General equations for the output from a nutating detector 
were known from Samuelsson's work, however these equations 
had not been solved for a specific detector. At the same 
time it was observed that the form of the covariance function 
for the noise offered an easy solution to integral equation 
of the Karhunen- Loeve expansion which made statistical 
detection theory appealing. Some work had been done on th 
form of temporal processor using statistical detection the >ry 
but this was limited to only one nutation because of the 
common background noise component between nutations. 

Harger's derivations provided a means for describing a 
temporal processor using statistical detection theory that 
based its decisions on multiple observations and included 
the background noise. His work was extended here to include 
unknown parameters, amplitude and position, which more 
closely characterized the system. This extension, called a 
threshold detector, was shown to be optimum in the case of 
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small signal-to-noise ratios where detection is most diffi- 
cult. The derivation also used the Gaussian assumption to 
describe the background. This may or may not be correct. 

But in practice, the detector's performance may not suffer 
greatly if this assumption is wrong. 

To implement the threshold detector, the integrals 
describing the detecor output were solved using the Gaussian 
quadrature method of numerical integration. Checking the 
computed coefficients by a summation provided a means of 
accepting the validity of the integrations. The form of the 
threshold detection system using a rectangular detector was 
determined and the frequency spectrum of the optimum filter 
was shown. To specify performance, the probability of 
detection is plotted against signal-to-noise and background- 
to-noise ratios and the number of nutations on which the 
decision was based. The signal spectrum for a point target 
and background power spectral density is also plotted. 

A designer may find these results useful in developing a 
system or investigating the performance of an actual system as 
compared to the optimum. He may also use these results to 
determine the size of a rectangular spatial filter to be used. 

Future work should include extending these results to 
spatial filters that are circular or elliptic. The same 
numerical integration algorithm could probably be used. The 
threshold detector was found to be highly sensitive to back- 
ground correlation lengths. Thus, some form of adaptive 
processor should be included to minimize the effect of 
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background. At the present time, a measurement program is 
underway at Naval Weapons Center, China Lake, California 
to determine average background correlation lengths but for 
different environments, these correlation lengths could be 
expected to vary considerably. 
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APPENDIX A 



DERIVATION OF SIGNAL SPECTRAL COEFFICIENTS 



The radiance on the detector is periodic because of the 
nutation therefore it can he expressed in a Fourier series, 



H ( t) = Z H e 
v J n 

n 



inw t 
- o 



where 



T -inu) t 

J c\ 



F n =i J H (t) e 0 dt 



(1A) 



(2A) 



But , 

H (t) = Jn ' (r) t (r-p (t) ) d 2 r . (3A) 

Substituting (3A) in (2A) 

, 7 n fT -jnaj t 

H = dr N'(r) dt x(r-p(t) e (4A) 

n j - i 0 

The transmittance, t, may be expressed in terms of the 

inverse Fourier transform, assuming an infinite image plane, 

( , - j 2iTk- (r-p (t) ) 

t C r ‘P C t ) ) = Jd Z k x * (k) e ~ ~ • (5A) 

Rewriting the exponent, 

k * (r-p (t) ) = k * r - pk x costo Q t - pk y simo Q t C^A) 

Equation (4A) can be rewritten. 



88 



T 

•n = N ' Cr) y f dt /d 2 k T *(k)e 



- j 2 k • r 



j 2 ttAs in (oj t+ 0) -jnw t 
e e 0 



where 






A ' ,k x + k y 



0 = tan -1 k x /k y . 



Observe that 



N 



(k) = Jd 2 r N 1 (r) e 



j 2 7 r k • r 



thus 



f 2 ine 1 f T -jn(u t + 0) 

H = Jd z k N* Ck)T*(k)e jno ^ J dt e 0 

o 



• e 



jnAsin (o3 Q t + 0) 



Using the relation 



j u) = t / e ' jn1 ’ c - izsin * <*». 

1 o 



the expression for the coefficients becomes 



H = fd 2 k N ' (k) t* (k) e- ln0 J 2^p Jk 2 +k 2 
n J ~ ~ n \ x y 



Let 



then 



<t> = tan ky/kx 



j n e‘- in ' t ’ = e jn0 



Rewriting (11A) 



(7A) 



(8A) 



(9A) 



(10A) 



(HA) 
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H n = j n Jd 2 k N' (k) t* (k) e 



* jn<J> 



n 



2lTp 



J' 



k 2 + k 2 

x y 



(12A) 
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APPENDIX B 



DERIVATION OF BACKGROUND CORRELATION COEFFICIENTS 
h 

The n power spectrum component of the background at 
frequency mo Q is given by 



i r 

6 n'Ti R B^ e 

o 



jnoj 



dr 



(IB) 



But R-fx) is the correlation function of detector output 

D 

averaged with respect to the starting time. 



T 

R b (t) = I / E[B(t + i)B(t)]dt 



(2B) 



Using the expression for B(t) from Part II, 



R 



B 



(t) = y / dt E 

o 



Jn ' (r)x (r^-p (t+x) ) d^ 



’/n' (r 2 ) x(£ 2 -p (t)) 



, 2 

d r 2 



(3B) 



Interchanging the order of expectation and integration, 



B 



(t) = y f dt / d2 r i / d 2 r 2 T [rr£ (t+x) J 



• T [r2"ef t ^] E [ N ' ( ri )N ’ c r2 ) J 

Assuming stationary background, 

e[n* ( ri )N- (r 2 )] = 4>^(r 1 -r 2 ) 



(4B) 



(SB) 
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where is the background covariance function on the 

image plane. 

It is now convenient to express the integrand in terms 
of two-dimensional spatial Fourier transform. To begin with 

, j 2 Trk • r 7 

• T(r) = Jt( k)e ~ ~ d z k (6B) 



where t ( k) is the Fourier transform of the aperture function 

r(r) and k = (k , k ). Then 

~ ~ x y 

T [rre^ t+T ^] T Df2 



-/■ 



d 2 k 1 r 



00 



/ d 2 k 



2 T *(h )e 



j27T(ki-ri-k 2 -r 2 ) 



~ j 2 7T [jk x p (t + x) -k 2 p (x)J 



( 7B) 



Substituting (5B) and (7B) in (4B) and interchanging the 
order of integration, 

R b (t) = Jd 2 k 1 x (k x ) /d 2 k 2 x*(k 2 ) 

T -j2ir£k.. , p(t + T)-k 2 ' 

dt e 





2_ e j2,t ii'ri [j2 



/ d2 ri 



r 7 *i 27T l5?‘r? 

J d ll < ^B^l’^2- )e 



(8B) 



The last integral becomes, with R = r-^-r- 

- j 2irk ? * R" 



fd 2 R4>'(R) 



- j 2nk 2 • r^ 



(9B) 



But the bracket above is the Fourier transform of g(r) , or 
the Wiener spectrum of background on the image plane. Let 



W 



-j2irk 2 -R 2 

B ( I ; 2 } = d ? 



(10B) 
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If F (k) is the Fourier transform of the point spread func- 
tion and Wg(k) is the Wiener spectrum of background on the 
object plane, then 

Wg (k) = | F o Ck) | 2 W B (k) ( 1 1 B) 



This relation is analogous to the output noise spectrum 
expression due to a noise input to a linear system. In terms 
of Wg(k) the last line of (8B) is 

r j 2ir (k , -k ? ) • r, ? 

W«(k) Je ~ Z ~ i d Z r 1 = Ck)<5(k 1 -k 2 ) (12B) 

Substituting this into (8B) and performing the integration 
with respect to k ? , 

Rg (x ) = fd 2 k | x(k) | 2 W’(k) 

n r T j 2 irk * [ p C t+T ) -p (t) ] 

• jJ e ~ ~ ~ dt (13B) 



Since p(t) = (p cosw o t, p sinw o t) , 

k • [p (t+x) -p (t) ] = k x pcosw Q (t+T) + k y psinw Q (t+T) 

-k pcosw t - k psinw t 
x o y o 

= A(t)cosw t + B (t) s ino) Q t = c(t)cos (o> Q t - rp C t)) 



where 



A(t) = p[k x cosw o T + k y sinw o T - k x ] 

B (t) = p[-k x sinw o T + k y cosw o T - k y ] 

C (t) = ^/a 2 (t) + B 2 (t) = p ^2(1 - cosw 0 x) (k 2 + k Z ) 
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and 



>Kt) 



tan 



1 B (t) 

ATtT 



Therefore the last 
x /-T j 2 tt c 

T J e 

o 



integral of (13B) is 

COS (w T+^) 

dt = J 0 (2 ttc) 



( 1 4 B ) 



(15B) 



where J (•) is the zeroth order Bessel function. Therefore, 
o 

the correlation function R^(t) is given by 

r b ( t) = /l T(R) l 2 w b ( ^ J o( 4TrpR sin "T") d2 ^ (16B) 

where 

k = Jk 2 +k 2 and d 2 k = dk dk . 

— v x y ~x ~y 



Finally, the power spectrum coefficient B n is, from (IB) 



B n = /|t(k) | 2 W'(k)d 2 k 



. 1 



T 




sin 




- jnw x 
e 0 dr 



(17B) 



where a = 2vkp. The last integral becomes with x = (o q t/2 , 
ir 

J Q (2a sin x) cos2nx dx 




- j — J J (2a sin x)sin 2nn dx (18B) 

o 

The cosine integral equals J^(a) and the sine integral 
vanishes [2]. Therefore on using (18B) 
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(19B) 



6 n = /l T( ^^| 2 l F o^| 2 W B (k)j2(27T P k)d 2 k 



This is the general power spectrum expression for a nutating 

detector output. The power at frequency nto Q is expressed in 

terms of aperture transform r(k), optical transfer function 

F (k) , Wiener spectrum of background Wg(k) and nutation 

2 

radius through J n (27Tpk) 
k-plane . 



The integration is over the entire 
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APPENDIX C 



EXPECTATION OF A J2 



The essential calculation required to find 

A l 2 ( {z n } ) = E B [A 12 C(Z n }}B)] (1C) 

for Gaussian background B is the expectation 







T T ? 




1 = e b 


exp 


/ A(t) B (t) dt - c / B Z (t)dt 






0 0 


_ 



c = M/N 



and 



A(t) 



N 



M 

E 

n=l 









C2C) 



B(t) may be expanded in its Karhunen- Loeve (K-L) repre- 
sentation provided the mean square value 

E(B 2 (t)) < « (3C) 

and thus the integral equation 

T 

A.<f>.(t) = / R R (t,u) <t>.(u)du 0 <_ t £ T (4C) 

1 1 0 

can be solved. A^ are called the eigenvalues of the equation 
and <J>.(t) are called the eigenfunctions. 

Using Mercer's theorem, the correlation function of the 
background may be expressed as 
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oo 



( 5C) 



^ Vk (t ) ^^(u) 

k 

The K-L expansion of B(t) is 



B(t) = £ b k <t> k (t) 



with 



T 

b k = / B (t) 4>g C t) dt 0 < t £ T ( 6 C ) 

where the coefficients, bj, , are uncorrelated and because 
B(t) is Gaussian the coefficients are statistically inde- 
pendent with means y k and variances 
Also expanding A(t) , 



A(t) = £ a k <|) k (t) 
k 



wi th 



a k = / A(t) <t>{;(t)dt, 



0 < t < T. 



(7C) 



Rewriting in terms of the expansions 



I = E 



B 



exp 



[ CO CO CO oo 

> (i J l a n b k *1" a £ b P dt 

o n k n k 



J CO oo 



c / £ £ b n b£ dt 

o n k 



CSC) 



Us ing 



6 xy = f Ct) dt 
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I = E 



bn 



exp 



oo oo 

Z (-«* a b* + i a*b ) - c E b b* 
2 n n 2 n n' n n 

n n 



(9C) 



Using the independence of the coefficients 



I = n exp 



n 



a b* + 


i a*b - c 


b 


2 


2 n n 


2 n n 


n 





OO 

= n exp 
n 



a 






a 


2' 


1 n 1 


' E bn exp 


_ p 


b - 




4c 




n 2c 





(IOC) 



For Gaussian random variables 



E (wxx*) 



exp w m x m*/ (1 -2wA)J 



1 - 2wA 



where 



m x = E (x) 



A = E [ (x*™ x ) (x*-mj) 3 



(11C) 



and w is a constant. 

Let x = b - a /2c 

n n 



and 

then 



m = y - a /2c 
x n n 



E bn exp 



a 2 
n 



b - - 

n 2c 



exp {-c| y n -a n /2c| 2 (l + 2cA^J 



(1 + 2cA ) 
^ n - 



(12c) 



Substituting (12C) into (IOC), and separating terms 
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B 

I = It (1 + 2 cA ) 
n n 



exp 



1 co | a | a 

1 v 1 n 1 n 

2 1 + 2 cA 

n n 



exp 



od y*(y -a /c) 

c v n J 

2 Z. 1 + 2 cA 
n n 



exp 



c ^ y n (y n- a n /c ) 
2 n 1 + 2cX n 



( 13 C) 



This result can be written in closed form by using the 



definition A = (a,d> ) 
n ^ ’ Y n J 



1 : ' a n ' V 

2 L 1 + 2 cA 

n n 



/ dt,A(t ) / dt 2 A(t 2 )q(t 1 ,t 2 ) ( 14 C) 

o o 



where 



q(t 1 ,t 2 ) 



2M 

N 



£ 0 r 4 ) (t , ) <j>* (t- ) 

l + 2cA T n v l' Y n v - 2 J 
n n 



( 1 5 C ) 



is such that operating on q(t^,t 2 ) with 

f f VW + — : R N C T 1 > t 2 ) 

oo /M 



(•)dt 1 dt 7 ( 16 C) 



and using Mercer's theorem yields 



_N_ 

2 M 



q(t 1 ,t 2 ) + / R B (T 2 ,t 3 ) q(t 3 ,t 1 )dt 3 = RgCT-^t^ ( 17 C) 



Likewise , 



oo u * fij -a /c) + y (y*-a*/c) 
c £ M n u n n J n v n n / J 



n 



1 + 2 cA 



n 



\ S ( m b C t ) - ^i)h(t) dt 



( 18 C) 
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where 



h <o ' £ V‘) (19C > 

n n 

such that 

T 

Nh ( t) + 2M / R B (t 3 ,t) h(t 3 )dt 3 = Mm b (t) (20C) 

o 

Taking the logarithm of both sides of equation (13C) and 
using the equations above 



£n I = Z (1 + 2cA n ) 
n 



-1 



+ -rrr / dt, 

4M 1 

o 



I 



T 

/ dt, 
o 



M 



S' ( W - W> 

m= 1 



q ’ t 2^ 



f 

o 



N 



m b^ • M 



i ^ (Z t (t) - s t (t)) 



h (t) dt 



£n I = 2 (1 + 2bA n ) 

n 



T T 

+ 1 ii 



M M M M 

z z £ (ti) z z (t 2 ) - z s £ (t 1 ) z z m (t 2 ) 

£= 1 m=l 1=1 m=l 



M M M M 

^ ^ ^ S^Ctl) Z_ S m (t ? ) 



1 ni'" 2 j n_T Aj i 0—1 

m=l £ = 1 A- 1 



m= 1 



q(t 1 ,t 2 )dt 1 dt 2 



T 

/ m, (t ) h (t) dt 
o 

T M 

2 f 1 Z (Z (t) 
Mo £= 1 



S £ (t)) h ( t) dt 



(2 1C) 
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The likelihood ratio desired is 



A 12«V’ 



exp 




T 

/ 

o 



Z n (t)S n (t)dt 



l M T 2 

'tV S n^ dt 



I . 



(22C) 



Taking the logarithm of both sides and substituting the 
expression for £^(1), 



A l2 ({Z n }) = 



- M T 

£ Z / Z (t)S (t) dt 
N t n v J n v J 

n=l o 



1 M T 2 
--If S 4 ( t) d t 
N i n v ' 

1N n=l o 



♦ l * n (1 4 f V 1 



T T , M M 

♦ m 1 1 * z d'P z W 

o o £ = 1 m= 1 



M M M M 

z W \ W - \ W Z W 

l-l m=l m=l x=l 



M M x 

1 Z q ^1 ’ t 2^ ) dt l dt 2 

!.= 1 m=l 



T 7 T M 

-/ m b (t)h(t)dt + ^ / S (Z £ (t) - S £ (t))h(t)dt 
o o & - 1 



(23C) 
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APPENDIX D 



KARHUNEN -LOEVE REPRESENTATION OF EQUIVALENT 
“ SIGNAL -TO -NOISE RATIO 



The variance of the threshold statistic 



M r T 

Z J z n (t)q(t)dt 

n=l o 



W = 



where 



q(t) - ± 



f(t) - / 



f (u)q(u,t)du 



(ID) 



under hypotheses H q and is 



T 

VAR = —/ q 2 (t)dt + M 2 fdfq(t) fduq(u)R (t,u) (2D) 

z o y 

Substituting (ID) into the first term of the variance 
expression 



NM 

z 



f 1 q 2 (t)dt = — J' dt f 2 (t) - 2 f (t) / T f(u)q(t,u)du 



T T 

+ / / f(u)f (v)q(t ,u)q(t,v)dudv 
o o 



(3D) 



Using the Karhunen- Loeve representation 



NM 

z 



f < 2 Ct)dt - f 

o 



CO 00 f j 

Z Z f,f * / 

k I k «■ 



<J> k (t)<f>* (t)dt 



2M 



w _ AT 1 

2 l 2 f k f ! 777 IT— I *k ft) *S (t ) dt / * m Cu)^(u)du 
k Z 1 IT x n 0 0 
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oo oo oo oo 



+ 1111 f, f* 



2m , 2m , 

N A m N A 



v o m „ k l , 2 m , , , 2 m , 

k l m n 1 + — - A 1 + -77 - a o 

N m N n 



r T 

— - J 4 > m (t)<P*(t)dt 



Using 



T t 

/ 4 > k (u) 4 >*(u)dt / 4 > n (v) 4 >^(v)dt 



( 4 D) 



{ nk ' / 



then 



NM 

z 



T / 2 } 

J d 2 (t)dt • f £ V^l- 2 ^ 



2 M , 

N A k 



2M . 

N A k 



2M , \2n 

N A k 

1 + 2M a 

N k. 



= ZM v f f* 

N f r k k l , 2 M . 
k \1 + -rr A, 

N k. 



( 5 D) 



The second term may be expressed in a similar manner to 



give 



VAR = 



— Z f, f* 
N k k k 



l + f A V 

N k> 



♦ (f 



.2 00 _ _ 



1 



k fk % t ^ X 



-V 



Ak 



N k/ 



VAR ■ ¥ 2 V{ j 4 2 M ~ 
k 1 ~T A k 



(( J) 



The equivalent signal to noise ratio in Part II is 



d = 



M 

Z 

n=l 



i v i T 

E • a n / f(t)q(t)dt 



( 7 D) 



S VAR 



The K-L representation of the numerator is 



M . f T . _ M 

T. a I f (t)q(t)dt = Z a 
i n J i n 

n=l o n-I 



|z |f k | 2 /(l*f £ k ; 



. ( 8 D) 
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Z|ts> 



d may now be represented by dividing (8D) by the square root 
of (6D) 



d = -f 



M vii' |2< ri 2M-. ^ 

£ a n ^ £ kl /U *— V 
n=l 



— Z I f I 
N 1 k l 



2/ i + 2M 
' 1 N A k J 



d = 



M 

M * a n ^ 
n=l 



?ihi 2 /c i+ f V 



JWji 



(9D) 
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APPENDIX E 



NON- ZERO TERMS FOR NUMERICAL INTEGRATION 

Equation (5.5) for the radiance function for a point 
target may be written as 



I 



where 

a 

3 

c 

P 




/< 



2 2 
-x -y 



-savV 

e ' 



s inax s inBy 
x y 



•e-jn* 




dxdy 



w/ / 2a 
h / /Jo 



/2p/a 



/2 



a 






Q = — (r -y ) 

H a v y 3 o' 

and 

(p = tan y/x (IE) 

Let E(x,y) be that part of the integral that is even about 
both planar axis which is 



E(x,y) - £ T 

Using Euler's relation the integral becomes 
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I - fn (x ,y) [cos (Px+Qy ) + jsin(Px+Qy) ] tcosncf) - jsin nc}>]dxdy (3E) 



Using trigonometric identities 




i-1 1 



where 



= cos Px cos Qy cos n<J> 

Ty = -sinPx sin Qy cos n<j> 

Tj = -sin Px cos Qy sin n<j> 

= -cos Px sin Qy sin n<j> 

= -j sin Px cos Qy cos n<J> 

Tg = -j cos Px sin Qy cos n<j> 

Ty = -j cos Px cos Qy sin n<j> 

T8 = +j sin Px sin Qy sin n<}> (4E) 



To save computation time it is necessary to determine 
which of the eight integrals are always zero and for which n 
they are zero. 

Observe that 




n 



= R e (cos<{> + j sin <j>) 
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Likewise 



sin ncf> 




I 



m 



(x+jy) n 



(6E) 



The first term in brackets is always even and may be made 
part of E(x,y). The second term in brackets may be expanded 
by the binomial expansion with the general term 

(x) n ' k (jy) k . (7E) 



The first integral of (4E) becomes 

Ij = jE(x,y)cos Px cos Qy R e (x+jy) n dxdy (8E) 



For 1^ to be non- zero, Re(x+jy) n must be even in x and even 
in y. From (7E) this is true only of even n; for odd n, 1^ 
will always be zero. 

Examining the other seven integrals yields 



I 



r jE(x,y)[cos Px cos Qy cos n4> 

+ j sin Px sin Qy sin n<{>]dxdy 

jE(x,y)[cos Px sin Qy sin n<£> 

+ j sin Px cos Qy cos n<t>]dxdy 



n even 



n odd 
(9E) 



Observing that symmetry about the x- and y-axis is necessary 
for this integral to be non-zero and integrating only when 
the integral is non-zero permits one to integrate only over 
the first quadrant . 
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APPENDIX F 





SUMMATION CHECK 


The signal 


for a point target was shown in Part II to be 


H (to t) 
v o ' 


= J N' (r)r(r-p(t)) d 2 r . (IF) 



Assuming an infinite image plane, H(w Q t) may be written in 
terms of the Fourier transform as 



H (to t) 
^ o J 


( - j 2irk • p (t ) ? 

= J N' (k)t*(k)e ~ ~ d Z k (2 F) 



For a rectangular detector (2F) may be written as 



H(0) o t) 


r _-2« 2 o 2 (k 2 x *k 2 y l sinnkx sini[hky 

1 TT k TTk 

J x y 

r i2vk- (p (t) -r ) -j2ir(k x +k y ) ^ ( 3F ) 

• 1 e • e y 


Let 





x = /2 airkx 
y = /7 anky 



and noting the 


integral may be separated into twoparts 


H(w Q t) 


/2 

, 2 i— x(pcosw t - (r -x )) 

- 1 fe" x sin Kx e ° 0X0 dx 

71 J /2 a 

hy 

s 111 — /2 

1 r 2 /2a i — y(psinw t - (r -y 1) 

lf-y o y o' j 

• — /e 7 e 7 ay . 

TT J y 

(4F) 
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Since the integration is -°° to +“> only even integrands 
be non-zero. Rewriting (4F) 



H(w 0 t) 



= I f 

T J 



00 2 
e x 



s m 



wx 

— '(12. COS f— x[pCOSti> t 
X \0 O 



(W 1 . 



2 -f 

TT -F 



sm 



hy 



■y L /2 o f/2 r . . f . 

7 cos — y psmu t - fr -y j 

y \o o K y 1 o J 



The general form of the integral is 



2 . 



R 



- RCB) - f / e‘ u cos(cu) 

r\ y 



du . 



Differentiating with respect to B 



dR(B) 2 f°° -u 2 „ 0 j 

" dB ~~ ' = tt J e C0S ^ cos ^ U ' 



But the solution to (7F) is known [2]: 



dR(B) 

dB 



2/tt 



(B-C) 2 . (B+C) 2 



+ e 



Since the constant of integration R(0) = 0, 



R(B) = 



2/tt 



f 



(v-c) 



f 



. ( y+c ) 



dv + j e 
o 



dw 



Let 



P = 



v-c 



and 



Q = 



w+c 
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will 

dx 

])dy. 

(SF) 

(6F) 

(7F) 

(8F) 

(9F) 



then 



R(B) = 



B-C 

M 2 

/IT J- C/2 



2 ? 

L -pL 2 

e v dP + z 



v 7T J- 



B+C 

2 2 
e'Q dQ 



it J-C/2 



c r r B-C\ j. n _c fB+C^ 

Erf C-o- ) + Erf (-r— ) 



(1 OF) 



whe re 



? r z _ 2 

Erf(z) = / e ^ dy . 



" IT O 



(11F) 



Substituting the result (10F) into (5F) 



H( % t) - f 



Erf 



(f[ 



pcosw t - r +x 



-Erf 



/2 r 



1 + - 
°J .n 



/2a 



v/ 



pcosw t - r +x 
o x o 



}-—«) 

J /2a / 



+ Erf |— [psinw t - r +y 1 + — h 

at o y ' o J ^ 



Erf \2I [ps inw Q t • r y -y 0 



]' 



/2a , 



(12J ) 



The mean square voltage of the background may also be 
solved in a like manner. 

From Part V, the mean square voltage was shown to be 



~ e l = /l T ^)| 2 IV^I 2 



C 1 3 F ) 



which for a rectangular detector becomes 



. 2, , -4ir 2o2 k 2 

-2 „ f sln e x ,, 

e^. = a n -2a / ^ — * — ^ ^ dk 

B 6 J u 2 k 2 a 2 + (2itk ) 2 * 

X X 
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e 



(14F) 



•a B -2B 



/ 



sin (irkyh) 



- 2 2,2 

-4ir o k 

y 



2,2 

TT k 

y 



6 



” + ( 2 tt ky} 



dk 

2 y 



Letting 



t = 2iTk 



the general form of the integral to be solved is 



. 2 



R(p) = - J sin (nt/2) e 



2 2 
-at 



t2/4 + t^ 



Taking the first and second derivatives one has 



and 



dR (y ) 1_ f sin (2y t) e 

du tt J t 



2.2 
-a t 



2 , . 2 
a + 1 



dt 



d R (y ) 
2 



2.2 



du 



2_ f cos (2y t) -a^t 

tt -L 2 2 

0 a + t 



The solution to the second derivative is known [2] 

- sLL fe-w Erfc (ao- f-) 
du 2 2 “ L 20 



where 



e ay Erfc (a a + ] 



00 2 

Erfc(z) = — J e ^ dy 
/tt" z 



is the complimentary error function. 

With lengthy and tedious calculations it can be 
by using the relation 



(15F) 



(16F) 



( 1 7 F ) 



( 1 8 F ) 



(19F) 



shown 
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f v 


J dv 


J F(w)dw 


0 


_ 0 



(2 OF) 



with initial conditions equal to zero and 



F(w) = -- O il 
du Z 



(21F) 



that 



R(y) 



a 



1 - Erfc 2° 



/ttu 



1 - e 



(y/2o) 



2 2 
a a 



2a* 



e 01 ^ Erfc (aa + ^-) 



+ e ay Erfc ' (aa - ^-) 

- 2 Erfc (aa) 

The mean square voltage thus becomes 
e^ = 4aBOg R(w;a)R(h;3) . 



(22F) 



(231 
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APPENDIX G - CALCULATED COEFFICIENTS 
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^ fc)r-«o v jcA<t-LAvOr^coo , ' 0 '~‘fMcn> 4 "*A\'Or--oooo'— '^Jco^LAvor^ coo^Or-KMroNtLAsor^coo^ 
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OOUUCOUOOUCOOOUCOUOUOCUUUOCCCCCUCU 



APPENDIX H 



COMPUTER PROGRAM TO CALCULATE THE FOURIER COEFFICIENTS 
OF THE BACKGROUND CORRELATION FUNCTION, A POINT TARGET 
AND A\| AVERAGED TARGET OVER A GAUSSIAN POINTING ERROR 
FOR A NUTATING OPTICAL SYSTEM WITH A RECTANGULAR 
DETECTOR. 



AA = BACKGROUND CORRELATION FUNCTION COEFFIZERO ORDER) 
ALFA = BACKGROUND CORRELATION COEFFICIENTS (VECTOR) 

F NO = SIGNAL FOURIER COEFFICIENT (ZERO ORDER) 

FN = SIGNAL FOURIER COEFFICIENTS (VECTOR) 

FAVGO = AVERAGED SIGNAL COEFF (ZERO ORDER) 

FAVG = AVERAGED SIGNAL FOURIER COEFF (VECTOR) 

RHO = NUTATION RADIUS 

SIG = STANDARD DEVIATION OF BLUR CIRCLE 
WT = WIDTH OF DETECTOR ( X-D I RECT I ON) 

H = HEIGHT OF DETECTOR ( Y-D IR EC T ION ) 

AL = INVERSE CORR LENGTH OF BACKGROUND ( X- D IR ECT I ON ) 

BE = INVERSE CDRR LENGTH OF BACKGROUND ( Y-DI RECT I ON) 

PO = STANDARD DEVISTIOM OF POINTING ERROR 

XO = DFTrCTOR COORDINATE X-DIRECTION 

YO = DETECTOR COORDINATE Y-OIRECTIGN 

RX = POINT TARGET COORDINATE X-DIRECTION 

RY = POINT TARGET COORDINATE Y-DI RECT I ON 

R.RANGLE = POLAR COORDINATES OF POINT TARGET 

NS = NUMBER OF FIRST COEFFICIENT DESIRED 

NUM = NUMBER OF LAST COEFFICIENT DESIRED 

I AVE= 1 " S IGNAL" WILL CALCULATE AVERAGE SIGNAL COEFF 

I A VE = 0 "SIGNAL" WILL CALCULATE POINT TARGET COEFF 



IMPLICIT REAL-V8 (A-H.O-Z) 

DIMFNSICN A L F A ( 120) 

COMPLEX^: 16 FN ( 120 ), FN 0, FAVG ( 12 0 ) , FAVGO 

COMMON /DTPRM/ RHO , S I G, W I , H, AL , BE , RX , RY , XO , YO 

COMMON /AVPRM/ NUM.IAVE.PQ 

COMMON /START/ NS 

COMMON /ANGL/ RANGL 

RFAD 2 , RHO . SIG, W I , H.AL ,BE 

READ 2 , XO, YO. R, ANGL 

READ 3. NS, NUM, I A VE , PO 

RANGL = 0 .01 745D0 > AN GL 

RX = R*DCGS ( RANGL ) 

R Y=RA DS I M RANGL) 

RANGL = DAT AN? ( RY-YO. RX-XO) 

CALL SUMCK 

CALL BKGRN( ALFA , AA) 

CALL SIGNAL(FN.FNO) 

WRITE (9) RHO • S I G , W I , H » AL , BE » RX » RY » X0» YO , PO 
I A VE= I 

IF(IAVE.EO.l) RX= 0 . ODO 
IF(IAVE.EO.l) RY = 0 . 0 DO 
CALL SUMCK 

CALL S IGNAL( FAVG, FAVGO) 

WRITE (9) FMO, PN, AA, ALFA 
WRITE (9) FAVGO, FAVG 

1 FORMAT (1015) 

2 FORMA T ( 8D1 0. 0) 

3 FORMAT ( 31 5 , ID 10 . 0 ) 

4 FORMAT (8F10. 4 ) 

STOP 

END 
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SUBROUTINE BKGRN(VP,V) 



SUBROUTINE 8KGRN COMPUTES THE COEFFICIENTS OF THE 
BACKGROUND CORRELATION FUNCTION UP TO 120 ORDERS. 
HIGHER ORDERS MAY EASILY BE COMPUTED BY ONLY CHANGING 
THE DIMENSION STATEMENTS. 



ZP = VALUES OF RADII FOR DIFFERENT ANNULI 
Z = WORKING VECTOR FOR RADII 
MAXR = MAXIMUM RADIUS OF OUTER ANNULUS 
NIZ1 = DIMENSION OF ZP 

VP = WORKING VECTOR IN FIRST PART AND RETURN VECTOR 
IN THE SECOND PART 
V = ZERO ORDER COEFFICIENT 



IMPLICIT REAL -8 (A-H,0-Z) 

DIMENSION VP< 1) ,DN( 50, 120) 

DIMENSION DI ( 50) , ZP( 20) , Z( 50) 

DIMENSION R ( 24 ) , W(24) ,TH(25 ) 

COMMON / A V P R M / NUM,IAVE,PO 

COMMON / DT PRM/ R HO , S I G, W I , H , AL , B E , R X, R Y , X 0 , YO 
COMMON /RK PRM/ A,B,C,D 
COMMON /N FIRST/ NFIRST 
DATA MAXR/4/, NR/ 1/ , NI Zl/20/ 

DATA ZP/O.ODOiI .200,2. 400,5.500 , 8 . 6D0 , 1 1 . 8 DO , 15 . DO , 20 . 
IDO, 25.0 DO, 30. ODO, 40.0D0, 50. ODO, 60. ODO, 70. ODO, 80. ODO, 
290. ODO, 100.0 DO, 1 5 0 . 0 DO , 200 . 0 DO, 400. ODO/ 

N F I RST = 0 

PI = DARCOS ( - 1 . DO ) 



CALCULATE CONSTANTS FCR INTEGRAL 

CONST = ALTBE*SIGT*2/PI /PI 
A=WI/2. DO/S IG 
B =H/2 . DO/ S IG 
C = ( AL* : SIG)--'*: ; 2 
D = { BE J >S I G) **2 
E = RHO/SIG 



SET ANNULI FOP SPECIFIC PARAMETERS 

DO 150 I = 1 , N I Z 1 
150 Z(I)=ZP(I)/E 
I = 0 

I F ( Z( N I Z1 ) .GT .MAXR ) GO TO 210 
DZ = PI /NR/E 
MMM= 50- NI Z 1 
DO 200 1 = 1 ^ MM M 
ZU+NIZ1) = Z ( N I Z 1 ) + I T D Z 
I F ( Z( I +N I Z 1 ) . GT . MA XR ) GO TO 210 
200 CONTINUE 



MAX = NUMBER OF ANNULI 
210 MAX= I +N I Zl-1 



PRINT 4,RH0,SIG, H, WI , AL , BE , A , B , C , D , E , DZ 
PRINT 4, ( Z( I ) , 1= I, MAX ) 

DO 300 1=1, MAX 



COMPUTE RADI I , WEIGHTS AND ANGLES FOR EACH ANNULUS 
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CALL WAA(Z( I) ,Z( 1+1) ,R,TH,W,NPQ) 

C 

C 

N F I RST = 1 
DI(I) = 0 .DO 
DO 115 NN= 1 , NUM 
115 DN(I,NN) = 0. DO 
C 
C 

C COMPUTE FUNCTION FOR EACH RADIUS 

C 

DO 280 K=l , NPQ 

CALL 3ESSL ( R ( < ) A E , VP , VJN , NUM) 

RV = VJN’ ? >2 / DEX P ( P ( K )**2 ) 

E P = DE XP ( R ( K ) * * 2 ) 

DO 12 0 NN= 1 , NUM 
120 VP(NN) = VP( NN)**2/EP 
R EC = 0 . D 0 
C 

c 

C COMPUTE FUNCTION FOR EACH ANGLE 

C 

DO 270 J=1,NPQ 

270 REC = REC+FU(R(K)tDCOS(TH(J)) ,R(K)*DSIN(TH(J) )) 
V= 0 • DO 

DO 13 0 NN=1 , NUM 

13 0 DN( I , NN ) =DN( I ,NN) + REC< : VP( NNI W < K) 

280 OKI) = DI(I) + REC 5 VRV*W( K) 

300 CONTINUE 
C 
C 

C SUM VALUE FROM EACH ANNULUS 

C 

DO 400 1=1, MAX 
J = MAX +1-1 
4 00 V = V + D I { J ) 

V=V*CONST 

PRINT 3 , N ♦ V , D I ( M AX ) 

NBQ = 0 

PRINT 1 3, NBO , V 

PRINT 4, (DKI), 1=1, MAX) 

DO 500 NP= 1 , NUM 
VQM = 0 . DO 
VP( NP) =0. DO 
DO 600 1=1, MAX 
J= MAX +1-1 

VP(NP) =VP(NP) +DN(J ,NP) 

1 F ( DN( J,NP).GT. VQM) VQM=DN ( J , NP ) 

600 CONTINUE 

VP(NP) = VP( NP) A- C ON ST 
PRINT 13, NP, VP(NP ) 

PRINT 1 4 , NP , VQM, DN ( MAX , NP ) 

500 CONTINUE 

C 

C CALCULATE SUM FOR INTEGRATION CHECK 

C 

SUM=V 

DO 700 J = 1 , NUM 
I = N UM + 1 - J 

SUM=SUM+2.D0VVP( I ) 

700 continue 

PRINT 12, SUM 
RETURN 

1 FORMAT (1015) 

2 FORMAT ( 6E 10. 0) 

3 FORMAT ( 21 10, 1P2E 15.6) 

4 FORM/ T{ 1P8E15.6) 

12 FORMAT ( 5X, G20. 13 ) 

13 FORMAT (5X, 14, 6X, G20. 13 ) 

14 FORMA T ( 5X,I4,6X,G20.13 ,6X,G20. 13) 

16 FORMAT ( 3( 1PD20. 12) ) 

END 
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OOOOOOO 



FUNCT ION FU(X, Y ) 



FUNCTION FU IS AN AUXILIARY ROUTINE FOR BKGRN THAT 
COMPUTES THAT PART OF THE INTEGRAND NOT RADIALLY 
SYMMETRIC . 



IMPLICIT REALMS (A-H,0-Z) 

COMMON /BKPRM/ A,B,C,D 
P I = 1. DO 

IFIX.NE .O.DOJ P I = D S I N ( A"<i X ) / ( Al 1 X ) 
IF(Y.NE.O.DO) P 1 = P 1* DS I N ( B* Y ) / ( B-Y ) 
F U= P 1 T 2 / ( C + X'i‘ 2 ) /( C+Y**2) 

RETURN 

END 
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SUBROUTINE S I GNA L ( REC , RECO) 



SUBROUTINE SIGNAL COMPUTES THE FOURIER COEFFICIENTS OF 
A POINT TARGET OR AVERAGED COEFFICIENTS OVER A 
GAUSSIAN POINTING ERROR 



ZP = VAIUFS OF RADII FOR DIFFERENT ANNULI 

Z = WORKING VECTOR FOR RADII 

MAXR = MAXIMUM RADIUS OF OUTER ANNULUS 

NIZ1 = DIMENSION OF ZP 

REC = RETURN VECTOR OF COEFFICIENTS 

RECO - ZERO ORDER COEFFICIENT 



IMPLICIT R E A L $ 8 (A-H.O-Z) 

DIMENSION BN (120) 

D I ME NS I ON DI (50) , DN J ( 5 0 , 1 2 0 ) , DN ( 5 0 , 1 2 0 >,Z(50),ZP(20) 
DIMENSION R( 24) , W( 24) ,TH( 24) 

COMPLEX-- 16 REC( 1 ) » VEC( 120 ) , RECO, VS, F, WP 
COMMON /A VPRM/ NUM , I AVE , PO 

COMMON /DTPRM/ R HO , S I G, W I , H , AL , B E , RX , R Y , XO , YO 

COMMON /SIPRM/ A,B,C,PX,QX 

COMMON /NFIRST/ NFIRST 

COMMON /START/ NS 

COMMON /ANGL/ RANGL 

DATA ZP/0.0D0,1. 2D0,2. ADO, 5. 5C0 ,8.600,11 .800,15.00,20. 
IDO, 25 .000, 30.0D0, 40. ODO, 50. ODO, 60 . OD 0, 7 0 . OD 0 , 80. 00 0 , 
290.000, 100. ODO, 150. ODO, 20 0.0 DO, 400. 000/ 

DATA MA XR / 7/, NR / 1 / , N I Z 1 /20/ 

N F I PS T = 0 

VS=DC MP LX( 0. 00,-1 .DO) 

PI = DAPCO S ( - 1 . DO ) 

P 12 = PI/2. 00 
KKK= 0 



CALCULATE CONSTANTS FGR INTEGRAL 
G AMA= S I G 

IF(IAVE.eQ.l) G A MA =CS QRT ( S I GM 2 + PO** 2) 
SG = DS QRT ( 2 . DO ) *- GAMA 
A = W I / SG 
B=H/SG 

CONST = A ri B/ P I / P I 
C = l.DO 
PX=-XO 
QX =— Y 0 

IF(IAVE.EQ.O) PX=PX+RX 
IF(IAVE.EQ.O) 3 X = £5 X+R Y 
PX = PX t DS QRT ( 2 . DO ) / GA M A 
QX=QX^ D SORT (2 . DO) /GAMA 
E = 2.D04RH0/SG 



SET ANNULI FOR SPECIFIC PARAMETERS 

DZ = 0 . DO 
DO 150 I = 1 , N I Z 1 
150 Z(I)=ZP(I)/(2.D0*E) 

I = 0 

IF(Z(NIZ1) .GT.MAXRJGO TO 210 
CZ= 15 . DO" PI/E 
M M M = 5 0 - NI Z1 
DO 200 I = 1 , MMM 
Z ( I +N I Z 1 ) = Z(NIZ1) + I :k DZ 
I F ( Z( I + NI Zl) . GT. MAXR) GO TO 210 
200 CONTINUE 

PRINT 1 7, Z l I +N I Zl ) 
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MAX = NUMBER OF ANNULI 



210 MAX=I +N I Zl-1 
GO TO 212 

211 MA X = I 

212 PRINT 4 ,RHC , SIG ,WI ,H, AL ,BE , RX, RY ,X0, Y0, PX, QX 
215 PRINT 4, ( Z( I) , 1= 1, MAX ) 

220 DO 300 1=1, MAX 



COMPUTE RADI I , WEIGHTS AND ANGLES FOR EACH ANNULUS 

CALL WAA(Z( I),Z( I + 1 ) , R , TH, W , NPQ ) 

NF I RST = 1 
NPR=NPQ 
D 1 ( I ) = 0 .DO 
DO 115 NN = NS,NUM 
DN J ( I , N N ) = 0 . D 0 
115 DN Cl, NN ) = 0 . DO 

I F ( KKK. EO.l ) GO TO 300 
DC 280 K = 1 , N P Q 
RK2=R(K)**2 

IF(RK2.GT.174J GO TO 285 
E P = DEX P ( RK 2 ) 

IF(EP.GT.1.D 50) GO TO 285 
CALL B E SSL ( R ( K ) 4 E , BN , V JN , N UM) 

RV=VJN/ EP 
RECG=F ( R( K) , TH, 0 ) 



COMPUTE FUNCTION FOR EACH RADIUS 

REC0=REC0*RV1 W ( K ) 

PP=RECO 

D I ( I) =D I ( I )+PP 
DO 130 NN = NS , N'JM 
BN(N'N) =BN(NN) / EP 
REC(NN) = F( R< K) , Trl, NN) 

REC(NN) =REC(NN)*BN (NN )*W(K ) 

P P = REC ( NN) 

DN ( I , NN ) =DN( I , NN ) + PP 
PP = REC ( NN )•-. VS 
DN J ( I , N N ) = DN J ( I ,NN)+PP 
130 CONTINUE 
280 CONTINUE 
300 CONTINUE 



SUM VALUE FROM EACH ANNULUS 

Ml = 0 . DO 
VI = 0 . DO 
V 2 = 0 . DO 
DO 400 1=1, MAX 
J = MAX +1-1 

IF(DABS(DI(J)).GT.CABS(VZ)) VZ=DI(J) 
I F ( DI ( J ) . L E . 0 . DO ) VI = Vl + DKJ) 

400 I F ( D I ( J) .GT.O.DO) V 2 = V2 + DKJ) 

RE C0= DC MP L X { V 1 + V2 , O.DO) 

R£CO=RECO-v CONST 
N B0 = 0 

PRINT 4, ( D! ( I ) , 1= 1,MAX ) 

PRINT 19,NB0,REC0 
PRINT 14,NBQ,VZ,DI (MAX) 

DO 500 NP=NS , NJM 
VQM = O.DO 
VI = O.CO 
VI J=0 .DO 
V 2 = O.DO 
V 2J = 0 . D 0 
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DO 600 1=1, MAX 

J =MAX+ 1 -I 

1F(DN(J,NP).LT«0.D0) VI = Vl + DM J , NP) 

1 F ( DN J ( J, NP) .LT .0 .DO ) V1J=V 1J + DNJ ( J, NP ) 

I F( DN( J , NP) .GE. 0. DO) V2 =V2+ DN( J , NP ) 

I F ( DM J ( J,ND) .GE.O.DO) V 2J = V 2J + DN J ( J , NP ) 

IF(DABS(DN(J,NP) ) . GT . DABS ( V Q,M ) ) VQM=DN( J , NP ) 

600 CONTINUE 

P.EC(NP) =DCMPLX( V 1+V2, V1J + V2J ) 

REC(NP) =F.EC(NP) *CONST 
PRINT 1 9, N° , REC ( NP ) 

PRINT 14, NP, VQM, DN ( MAX, NP ) 

500 CONTINUE 
725 CONTINUE 

VS = DCMP LX ( 0 . DO , 1 . DO ) 

DO 800 NP = NS , NUM 
800 REC (NP ) =R EC< NP) *VS**NP 
1 F ( IAVE.EQ.0 ) GO TO 960 
SUM=RECO 

IF(NS.NE.l) SUM = 0 . ODO 
PRINT 1 9 , NBQ, RECO 
C 
C 

C CALCULATE SUM FOR INTEGRATION CHECK 

C 

DO 900 I = NS , NUM 
PQ=RE C ( I ) 

PRINT 1 9, I, REC(I) 

900 SUM=SUM+2.D0#PQ 

PRINT 1 8 , SUM , SUM J 
RETURN 

960 CONTINUE 

VS=DCMPLX( O.DO,-l.DO) 

PRINT 1 9 , NPQ , RECO 
SUM=RECO 

IF(NS.NE.l) SUM=O.DO 
SUMJ=0 . DO 
DO 950 NP = NS , NUM 

V EC ( NP ) =R EC (NP)*DCMPLX(DCOS(NP RANGE ) , D S I N ( NP*RA NG L ) ) 
PRINT 19, NP, VEC(NP) 

PP = VEC ( NP ) 

SUM=SUM+2 .DO*PP 
P P=VE C ( NP)*VS 
SUMJ = SUMJ +PP 
950 CONTINUE 

PRINT 18, SUM, SUMJ 
RETURN 
285 KKK= 1 

GO TO 300 

1 FORMAT (1015) 

2 FORMAT ( 8E 1 0. 0) 

3 FORMAT ( 21 10, 1P2E 15 .6) 

4 FORMAT ( 1P8E15.S) 

5 FORMAT (215) 

6 FORMAT ( 1P2E15 .6 ) 

12 FORMAT ( 5X,G20. 13) 

13 FORMAT ( 5X, 14, 6X,G20.13) 

14 FORMAT ( 5X , 1 4 , 6X , G20 . 1 3 , 6X , G 20 . 1 3 ) 

15 FORMA T( 2X,I3,2X,1PE15.6,2X,1PE15.6,2X,I3,2X,1PE15.6,2X 
1, 1PE15.6) 

16 FORMAT ( 3 ( 1P020. 12 ) ) 

17 FORMA T ( 5X , 'MAX NOT LARGE ENOUGH ZM AX = ' , FI 5 . 8 ) 

18 FORMAT ( 2X, 'SUM= ' , G20. 13, ' SUMJ = ' , G20.13) 

19 F ORMAT ( 5 X , I 4 , 6X , G2 5 . 1 3 , G20 . 1 3 ) 

END 
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FUNCTION F ( R » TH » NN ) 



FUNCTION F IS AN AUXILIARY ROUTINE C 0R SIGNAL THAT 
COMPUTES THAT PART OF THE INTEGRAND NOT RADIALLY 
SYMMETRIC 



IMPLICIT REAL-53 (A-H,0-Z) 

COMPLEX*- 16 F 
DIMENS ION T H ( 1) 

DI ME NS I ON AAE (2 4) , BBE(24) , A AO (24 ) , BB0I24 ) 
COMMON /SIPRM/ A,B,C,PX,QX 
DATA NPP/24/ 

JJ=(-l)*tNN 
A A= 0. 00 0 
BB=0 . 0 DO 

IF(NN.GT.O) GO TO 5 
DO 1 J=1,NPR 
X=R4DC0S(TH< J ) ) 

Y=R^D SI N( TH ( J ) ) 

DSX = DS I N( PX'-X ) 

DSY=DSI N( QXL Y) 

DCX=DCO S ( P X & X ) 

CC.Y = DCO S ( QX~' Y ) 

P 1=1. DO 

IF(X.NE.O.DO) P 1 = DS IN( AX X) / ( A*'X ) 
IF(Y.NFrO.DO) PI = PIX'DS I N ( B#Y ) / (8>Y) 

AAE( J ) = DC X >DC Y» P 1 
BBE { J ) = DS X-DS Y A P 1 
A AO ( J) =-DCX DS Y PI 
BBO( J ) = -DSX s >DCY : ’P 1 
1 CONTINUE 

5 IF(JJ.LT. 0) GO TO 10 

DO 20 J = 1 1 NP R 

AA=AA + AAE( J KDCOS (NN-*TH(J ) ) 

BB = BB + BBE( J)*DSI N ( NN* TH ( J ) ) 

20 CONTINUE 

F = DCM PLX ( AA t BB) 

IF(NN.EQ.O) F=DCMPLX( AAtO.DO) 

RETURN 

10 DO 40 J = 1 » NPR 

A A= AA +A AC ( J ) * D SI N ( NN* TH( JJJ 
BB=BB+BBO( J )x DCCS (NN-THIJ ) ) 

40 CONTINUE 

F=DCMPL X( AA, BB) 

RETURN 

END 
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SUBROUTINE WAA(Ri,R2,R,TH,W.N) 



SUBROUTINE WAA COMPUTFS THE RADIAL PO I NT S . WE IGHT S AND 
ANGLES AT WHICH THE INTEGRAND MUST BE EVALUATED. 



NP = DEGREE OF GAUSSIAN LA 
R = RADIAL POINTS TO BE USED 
TH = ANGLES TO BE USED 
W = WEIGHTS TO BE USED 

NFIRST PREVENTS COMPUTING G-L POINTS MORE THAN ONCt 



IMPLICIT REALR8 ( A-H.O-Z) 

DIMENSION R ( 1 ) » W ( 1 ) » TH ( 1 ) 

R FAL *6 RC( 24J , WC ( 24) 

CGMMQN /NFIRST/ NFIRST 
DATA N P / 2 4 / 

P I = OARC OS ( — I .DO) 

PI N=PI / 2 . CO/NP 
I FINFIRST.. NF. 0) GO TO 99 
CALL GL024 (RC. WC. 0 . DO. 1 .DO ) 

DO 10 I =1 . NP 

10 T H ( I ) = ( I-.5D0H-P IN 

NFIRST = 1 

99 CONTINUE 
N = NP 

T = R2 - R1 
TW = T*r ( R2 + R1 ) 

F = TW+PI/NP 
DO 100 1=1 .NP 

RX = R 1 
RX = RX^Rl 
RY = R2 
RY = R Y'.'R? 

R(I) = DSORTIRX + RC(I)*(RY - RX)) 
W ( I ) = WC { I ) ?» F 

100 CONTINUE 
RETURN 
END 



124 



*w -H •£- *•> ^ 



SUBROUTINE GL 024( X , A, C , D ) 

C 

c 

C GAUSSI AN-LEGENDER INTEGRATION ROUTINE (24 POINTS) 

C 

r. 

DOUBLE PRECISION CMC, DPC 

DOUBLE PRECISION C,D,X( 24), A( 24),XX( 12), AA ( 12) 

DATA XX ( 1 ) , A A ( i ) , XX (2 ) . AA. ( 2 ) . XX ( 3) , AA( 3 ) , XX(4) , AA( 4) ,X 

1 XX< 6) , AA ( 6) ,XX( 7) , AA ( 7 ) . XX ( 8 ) , A A ( 8 ) , XX ( 9 ) , AA (9 ) , XX ( 10 

2 XX( 11 ) , AA( 11) , XX( 12) , AA( 12) / 

*.9951 87 2 19997 02 1360 17999740 97 D-0, . 1234122 9 7 99987199540 
* .974728 555 971 30949 8 19 8391993 OD-O ,. 2 853 1 3 88628933663 181 
v .93827455200273275852364900170-0, .44277438817419806163 
4-. 8864 15 52 7004401 0342 1 3 1 5434 1 9 D-0 , .59 29 8 5 849 15 4367 3 0 746 
* .82000 1985 9739029 2 1953949 8726D-0, . 73346481411080305739 
. 7401 24 1915 78554 3642438 28 10 31C-0, . 86 19 0 1 6 1 5 3 1 9532 7 59 1 7 
. 648093 651 93 697 5 56 92 52495 78o9 D-0 , . 976 1 8 6 5 2 104 1 138 8 82 69 
. 54 542 14 71 38883 95 3 5658 37561 72D-0,. 10744427011596563473 
.433793507626045 136 467 0842 3 19D-0, .11550566805372560135 
. 31 5042. 679696 163 3 7 4386 79329 13D-0, . 12 1 6 7 0472 9278033 9 1 20 
.191) 1 8 86747 36 16 309 15 8639 82 07D-0, .12583745634682829612 
. 64056 8 92 8626056 26 08 5043 082 62 D- 1 , .127938 19534675215697 
DMC = . 5D0-UD-C) 

DPC = . 5D0*(D+C) 

DO 2 1=1,12 
NI = 25 - I 

X(I) = —DMC* XX ( I ) + DPC 
X ( N I ) = DMC- XX(I) + DPC 

A (I ) = DMC--AA( I ) 

2 A ( N I ) = DMC- AA(I) 

RFTURN 
END 
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SUBROUTINE RESSL ( AG, P.F ,BZ , NUM) 



SlJR ROUT INF BFSSL COMPUTES BESSEL FUNCTIONS CF THE 
RECEIVED ARGUMENT 



AG - ARGUMENT OF BESSEL FUNCTION 

BF = RETURN VECTOR OF BESSEL FUNCTIONS 

BZ = ZERO ORDER BESSEL FUNCTION 

NIJM = HIGHEST ORDER OF BESSEL FUNCTION TO BE RETURNED 
•NY. NX = ORDFRS CALCULATED BY ASYMPTOTIC EXPANSION 



IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION BF3( 1000) ,BF ( 1) 

NZ = 1000 
KT=0 

IFIAG.LT. 0) KT= 1 

IFIAG.LT . 0 . DO ) AG=DABS ( AG) 

X =AG 

N Y= I F I X ( SNGL ( AG ) ) 

IF(NY.LE.L) NF- L6 
I F ( NY . F 0. ?. ) NF = 2 0 
I F ( NY . EQ. 3)NF=24 

IF( NY. GE. 4. AND.NY.LT. 6) N F= 6E-N Y 
I F ( NY . G F . 6 . AND . NY.LT.10) NF = 5P NY 
IF( NY. GE.10.AND.NY.LT. 15) NF = 44 NY 
I F ( NY . G E . 1 5 . AND . NY . LT . 29 ) N F = 3* N Y 
IF(NY.GE.29.AND. NY. LT. 80) NF = 2v ! NY 

IF( NY. GF.80.AND.NY.LT. 150) NF=IFIX(1.5^FLOAT(NY)) 

I F < NY . G E . 1 50 . AND . N Y . L T . 33 0 ) N F= I F I X ( 1 . 4 ; v FLOAT ( NY ) ) 
IFINY.GE. 330) NF = I F I X( 1 . 25'TF LOA T ( NY) ) 

I F< I FIX ( FLOAT (NF ) / 2 . ) . L T . I F I X ( F LG AT ( NF-»- 1 ) / 2 . ) ) NF=NF+ 1 
IF(NY.LT.L) NF = 8 
NY=NF 
NX=NY- 1 

SFT ORDERS > NY TO ZERO 

DO 1 I = NY . N7 
BF3( I ) = 0. D 0 
1 CONTINUE 



COMPUTE ASYMPTOTIC EXPANSION 

DO 2 MlJ =N X » NY 
Z=X/ DFL OAT (MU) 

AMU=DFLGAT (MU) 

A = DSORT ( l.DO-Zft-2) 

A1 =1 .500* (DLuG( ( l.DO+AJ/ZJ-A) 

DEL = DEXP( 2. DO/3 . DOpDLOG ( A1 ) ) 

XI = DEXP( 2.D0/3.D0^DLGG(4MU)+DLGG(DEL) ) 

CALL TABLE ( X 1 , A I V . A I PV ) 

D 1 = DEXP( 1 . DO / 3 . D O'*- DLOG < A MU ) ) 

D2 = OEXP ( 5 . DO / 3 . DO--DLOG( AMU) ) 

BO = -5 . DO/ (43 . D0? DEL»- + 2 ) + ( L .DO/ DSORT ( DEL ) )-/( 5.00/ ( 24. 
1D0 : ? A 4? 3 )- 1 . 0D0/( 3. ODO*- A ) ) 

FI = DtXP(0.25D0*DL0G<4.D0*-DcL ) -0 . 5D0>D LOG ( A ) ) 

B F 3 ( M U ) =F I M AI V/DL+AIPV-B0/D2 ) 

? CONTINUE 
NY2=NY-2 



USE RECURSIVE EQUATION TO GENERATE ORDERS < NX 

DO 4 1=1, NY? 

N = NY 2— I +1 

BF3 ( N) = (2.D0-DFL0AT (N+i ) /X )->B F3 ( N + l )-BF3(N+2) 
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4 CONTINUE 

SUM = 2 . D0 : ?BF3 I 1 ) / X— 8 F 3(2) 

DC 76 I = 2 • NY2 » 2 

SUM = SUM 6 2 , DO z 8 F3 ( I ) 

76 CONTINUE 



COMPUTE NORMALIZING CONSTANT 

SNORM = 1. DO/SUM 

DO 77 1=1, NUM 

ft F ( I ) =S NOR Mf BF3 (I) 



OCD ORDER BESSEL FUNCTIONS ARE ODD 

IFIKT .EO.l) BF( I)=BF( I )* ( - 1 .DO) $+ I 
77 CONTINUE 

B7. = ( 2.D0*RF3( 1) /X-BF3<2) )*SNORM 

RETURN 

END 
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SUBROUTINE TABLE ( X , AI V, AI PV) 



SUBROUTINE TABLE IS AIM AUXILIARY ROUTINE FOR BESSL 
THAT COMPUTES AIRY FUNCTIONS 



IMPLICIT REAL48 (A-H.O-ZJ 

DIMENSION A I (50 ) , A I P ( 50 ), ABI ( 15 I , ABIP( 15) 

DATA A I / 4*0. ODO, .56228000,4*0.000, . 560462 DO , 4-- : 0 . 0 DO , 

1 . 5 5872 4D0, 4' 0. OD 0, . 5 57058D 0 , 4* 0 . ODO , . 5 5545 6D0 ,44 0 . ODO , 
2. 553912 DO .4": 0 .ODO, . 55242 IDO , 4*0 . ODO, . 55098000, 4*0. ODO, 
3. 549504D0,4r 0. ODO , . 548 2 30 DO/ 

DATA AI P/4-.0. ODO, . 5668 73D 0, 4* 0. ODO ,. 56944 8D 0 , A - * 0 . ODO , 
1. 57192700,450. ODO, .574320 00,440 .0 DO, .57663500, 4*0 . ODO. 
2. 578878D0, 4-0. ODO, . 58 1 056D0 , 4* 0 . ODO , . 533 1 74D0 , 4* 0 . 0 DO , 
3 . 585235 DO . 4- 0 .ODO, . 58724500/ 

DATA AB 1/430. ODO, .5482 3 0 DO, . 5456 36 DO , . 543 180 DO , 

1 . 540844D0, .53861800, . 53648SD0,. 534448D0, . 532488D0, 
2.5306 01 DO, .5287 8 3 00, .52702700/ 

DATA A 8 ] P/440. ODO. . 587245 DO , . 59 1 12 ODO, .594823 DO, 

1 . 59837 2 DO, .60 178 2D 0, . 60 506 8D0, . 608239D0. . 61130500 , 

2. 6 142 75 DO, .6 1715600, .6199 54 00/ 

L = 0 

Y - DEXP ( DLDG( 1. 5D0 )-l . 5D0*DLDG( X) ) 

W = DEXP(DLOG(2.DO/3.DO)+1.5DO*DLOG(X)) 

L= 1 

I F ( Y . GT .1 . 5D0 ) GO TO 14 
I FtY.GT.O. 49999 9 DO ) GO TO 1 
L = 2 

W1 = 10 .DO * Y 
I Wl = I FI X( SNGL ( W1 ) ) 

W2=W1+1 . 

I W2 = I FIX! SNGL (W2 I J 
I 7 = I W 1 3 1 0 

I F C IW1.LT.IW2) I Z = I Z+ 5 
I Y = I Z - 5 

I E( I Y. FO. 0) GO TO 13 
GO TO 8 

1 Y = 10. DO Y 

I Y= I F I X ( SNGL ( Y) ) 

IF(IY-Y) 2,4,6 

2 I 7 = I Y + 1 
GO TO 8 

4 AIV=AI ( IY ) 

A I PV= A I P( I Y) 

GO TO 19 
6 I Z = I Y 

I Y= I Y- 1 

8 DP=Y-FLGAT( IY) 

IFIL.E0.2) DP= DP / 5 . DO 
IF(L.EQ.l) GO TO 12 
A I V = ( 1. DO-DP ) t A I ( I Y I+DP-sA I( I Z ) 

A I P V = (l.DO — DP)*AIP( IY)+DP*AIP< I Z J 
GO TO 19 

12 A I V = ( l.DO-DP)«-ABI(IY)+DP*ABI<IZ) 

A I P V = (1, . DO-OP )*ABIP( IYl+DP-'-AB IP( IZ) 

GO TO 1 9 

13 DP = Y/5.D0 

A I V = ( 1 . DO -DP )« C 0.564190D0 + DP-‘AI( IZ) 

A I PV = ( 1. DO-DP) *0.564190D0 + DP**AIP( IZ) 

19 A I V = DEXP(-DLOG ( 2 .DO )-0 . 2 5 D0*D LOG ( X )-W+DLOG( AI V) ) 

A I PV = DrXP(0.25D0- ; DL0G(X)-W + DL0G(AIPV) ) 

A I P V = A I P V (- .500) 

RETURN 

14 PRINT 30 

30 FORMAT ( 5X, 'HmVE EXCEEDED TABLE VALUES 1 ) 

RETURN 

END 
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SUBROUT INE SUMCK 



SUBROUTINE SUMCK COMPUTES THE ANALYTICAL EXPRESSION 
FOR THE INFINITE SUM CF COEFFICIENTS 

SISUM = INFINITE SUM OF POINT TARGET (IAVE=1) OR 

AVERAGED TARGET ( I AVE=0) COEFFICIENTS 

BKSUM = INFINITE SUM OF BACKGROUND COEFFICIENTS 



IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /DTPRK/ RHO , SI G , WI . H , AL , BE , RX , RY , XO , YO 
COMMON / A VPR M/ NUM.IAVE.PO 
COMMON /ANGL/ RANGL 
W = WI 

CONST =0.2500D0 

GAMA= DS ORT (SI G* R2+PO^*2) 

1F( IAVE.EO.O) GA MA = SI G 
DIV=2 .QO-OSQRT (2.DOKGAMA 
I F ( IAVE.EO.O) GO TO 10 



COMPUTE SUM OF AVERAGED COEFFICIENTS 

A=DERF( ( 2..D0 M XO + RHO )+W)/DI V) 

B^DER F ( ( 2 ..DO-' ( XO+RHO) -K )/ D IV ) 

C =DER F ( I 2 . 00 OF ( Y 0 ) +H ) / D I V) 

D= DER F ( (2.0n0*(Y0)-H) /DIV) 

4 SI SUM = CONST ■ I A- B ) $ (C-O) 

RI=R( W.AL.SIG) 

R2=R ( H . BE. SI G) 



COMPUTE SUM OF BACKGROUND COEFFICIENTS 

B KS UM = 4 . D 0 A A L I- 8 E R 1 R 2 / I W*H) F * 2 

PRINT 5.BKSUM.SISUM 

RETURN 

10 CONTINUE 



COMPUTE SUM OF POINT TARGET COEFFICIENTS 

A=DERF ( (2..0D0* (XO+RHO-DCOS I RANGL )-RX)+W ) /DI V) 
R=DERF ( I 2 . 00 OF { X0+ RHOF DCGS ( RANG L ) —RX ) -W )/ CIV) 
C= DER F ( ( 2 .D0£( YO+RHO* OSIN ( RANGL ) -RY) +H) /DIV) 
D=DERF ( (2.DO*(YO + RHO*'DS INI RANGL )-RY )-H) /DI V) 
GO TO 4 

5 FORMAT! 2X , 'BKSUM= ',015.8, 10 X, *SISUM=' .015.8) 
END 
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FUNCTION R ( II . Z . S I G ) 



FUNCTION R IS AN AUXILIARY ROUTINE FOR BKSUM 



IMPLICIT RE£L*8 (A-H.O-Z) 

P I =DA RC OS ( — 1 .DO) 

R=0 .DO 

F 1 =- ( ( U / ( 2 . DOF-S I G ) ) * + 2 ) 

F 1 = DFXP (FI) 

F 2 = ( Z S I G ) • 2 

E 2=DEXP ( F2 ) 

F 31=2*0 

E3 = DEX P (FBI ) 

F4=DE X P ( — E3 ) ) 

X1 = IJ/ ( 2.D0. SIG) 

YL =1 . DO-DERF (XI ) 

X2 = Z ; F S I G+ U / ( 2. DOFSIG) 

Y 2= 1 . DO-DFRF( X2) 

X3-Z1-SI G-U/(2. DOiS IG) 

Y3= 1 . D O-DE RF ( X3) 

X4=Z*S IG 

Y4=l. DO-DERF(XA) 

R = ( 0 / ( Z')v2) i ' ( 1. DO-Yl-( SI 3 w 2. 0 0/(0 SORT ( PI i*U) )* (1 . 00-E 
R=R+( E2/ (2 .D0*Z*-*3J )*< E3* Y 2 + E 4<- Y 3- 2. DO * Y A ) 

RETURN 

END 
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